
<!DOCTYPE html>
<html lang="en">
  <head>
    <meta charset="utf-8"> 
    <meta name="viewport" content="width=device-width, initial-scale=1">




<link
  href="http://maxcdn.bootstrapcdn.com/bootstrap/3.3.4/css/bootstrap.min.css"
  rel="stylesheet">

<style>
.julia {font-family: "Source Code Pro";
        color:#0033CC;
        }
body { padding-top: 60px; }
h5:before {content:"\2746\ ";}
h6:before {content:"\2742\ ";}
pre {display: block;}
</style>

<script src="http://code.jquery.com/jquery.js"></script>
<script src="http://maxcdn.bootstrapcdn.com/bootstrap/3.3.4/js/bootstrap.min.js"></script>

<!-- .julia:before {content: "julia> "} -->

<style></style>

<script src="https://cdn.plot.ly/plotly-latest.min.js"></script>


<!-- not TeX-AMS-MML_HTMLorMML-->
<script type="text/javascript"
  src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_SVG">  
</script>
<script>
MathJax.Hub.Config({
  tex2jax: {
    inlineMath: [ ["\$","\$"], ["\\(","\\)"]]
  },
  displayAlign: "left",
  displayIndent: "5%"
});
</script>


<script type="text/javascript">
$( document ).ready(function() {
  $("h1").each(function(index) { 
       var title = $( this ).text()
       $("#page_title").html("<strong>" + title + "</strong>");
       document.title = title
  });
  $( "h2" ).each(function( index ) {
    var nm =  $( this ).text();                                    
    var id = $.trim(nm).replace(/ /g,'');
    this.id = id
    $("#page_dropdown").append("<li><a href='#" + id + "'>" + nm + "</a></li>");
  });
  $('[data-toggle="popover"]').popover();  
});
</script>

</head>


<body data-spy="scroll" >

<nav class="navbar navbar-default  navbar-fixed-top">
  <div class="container-fluid">
    <!-- Brand and toggle get grouped for better mobile display -->
    <div class="navbar-header">
      <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#bs-example-navbar-collapse-1">
        <span class="sr-only">Toggle navigation</span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
      </button>
         
    </div>

    <!-- Collect the nav links, forms, and other content for toggling -->
    <div class="collapse navbar-collapse" id="bs-example-navbar-collapse-1">
      <ul class="nav navbar-nav">
        <li><a href="#" id="page_title"></a></li>
      </ul>
      <ul class="nav navbar-nav navbar-right">
         <li class="dropdown">
           <a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
           Jump to... <span class="caret"></span></a>
          <ul class="dropdown-menu" role="menu" id="page_dropdown"></ul>
        </li>
      </ul>
    </div><!-- /.navbar-collapse -->
  </div><!-- /.container-fluid -->
</nav>

<header>
</header>

<div class="container-fluid">
  <div class="span10 offset1">
<h1>A SymPy tutorial</h1><p>This tutorial provides an introduction to using <code>SymPy</code> within <code>Julia</code>. It owes an enormous debt to the tutorial for using SymPy within Python which may be found <a href="http://docs.sympy.org/dev/tutorial/index.html">here</a>. The overall structure and many examples are taken from that work, with adjustments and additions to illustrate the differences due to using <code>SymPy</code> within <code>Julia</code>.</p><p>This tutorial can be read as an <code>IJulia</code> notebook <a href="http://nbviewer.ipython.org/github/jverzani/SymPy.jl/blob/master/examples/tutorial.ipynb">here</a>.</p><p>After installing <code>SymPy</code>, which is discussed in the package's <code>README</code> file, we must first load it into <code>Julia</code> with the standard command <code>using</code>:</p><pre class="sourceCode julia">using SymPy</pre>
<p>The start up time is a bit lengthy.</p><h2>Symbols</h2><p>At the core of <code>SymPy</code> is the introduction of symbolic variables that differ quite a bit from <code>Julia</code>'s variables. Symbolic variables do not immediately evaluate to a value, rather the "symbolicness" propagates when interacted with. To keep things manageable, SymPy does some simplifications along the way.</p><p>Symbolic expressions are primarily of the <code>Sym</code> type and can be constructed in the standard way:</p><pre class="sourceCode julia">x = Sym("x")</pre>
<div class="well well-sm">
$$x$$</div>

<p>This creates a symbolic object <code>x</code>, which can be manipulated through further function calls.</p><p>There is the <code>&#64;syms</code> macro that makes creating multiple variables a bit less typing, as it creates variables in the local scope – no assignment is necessary. Compare these similar ways to create symbolic variables:</p><pre class="sourceCode julia">@syms a b c
a,b,c = Sym("a,b,c")</pre>
<pre class="output">
(a,b,c)</pre>

<h3>Assumptions</h3><p>Finally, there is the <code>symbols</code> constructor for producing symbolic objects. With <code>symbols</code> it is possible to pass assumptions onto the variables. A list of possible assumptions is <a href="http://docs.sympy.org/dev/modules/core.html#module-sympy.core.assumptions">here</a>. Some examples are:</p><pre class="sourceCode julia">u = symbols("u")
x = symbols("x", real=true)
y1, y2 = symbols("y1, y2", positive=true)
alpha = symbols("alpha", integer=true, positive=true)</pre>
<div class="well well-sm">
$$\alpha$$</div>

<p>As seen, the <code>symbols</code> function can be used to make one or more variables with zero, one or more assumptions. </p><p>We jump ahead for a second to illustrate, but here we see that <code>solve</code> will respect these assumptions, by failing to find solutions to these equations:</p><pre class="sourceCode julia">solve(x^2 + 1)   # ±i are not real</pre>
<div class="well well-sm">
\begin{bmatrix}\end{bmatrix}</div>

<pre class="sourceCode julia">solve(y1 + 1)    # -1 is not positive</pre>
<div class="well well-sm">
\begin{bmatrix}\end{bmatrix}</div>

<p>The <code>&#64;syms</code> macro can also have assumptions passed in as follows:</p><pre class="sourceCode julia">@syms u1 positive=true u2 positive=true
solve(u1 + u2)  # empty, though solving u1 - u2 is not.</pre>
<pre class="output">
0-element Array{Dict{SymPy.Sym,SymPy.Sym},1}</pre>

<p>As can be seen, there are several ways to create symbolic values. One caveat is that one can't use <code>Sym</code> to create a variable from a function name in Base.</p><h3>Special constants</h3><p><code>Julia</code> has its math constants, like <code>pi</code> and <code>e</code>, <code>SymPy</code> as well. A few of these have <code>Julia</code> counterparts provided by <code>SymPy</code>. For example, these three constants are defined (where <code>oo</code> is for infinity):</p><pre class="sourceCode julia">PI, E, oo</pre>
<pre class="output">
(pi,E,oo)</pre>

<p>(The pretty printing of SymPy objects does not work for tuples.)</p><p>Numeric values themselves can be symbolic. This example shows the difference. The first <code>asin</code> call dispatches to <code>Julia</code>'s <code>asin</code> function, the second to <code>SymPy</code>'s:</p><pre class="sourceCode julia">[asin(1), asin(Sym(1))]</pre>
<div class="well well-sm">
\begin{bmatrix}1.5707963267949\\\frac{\pi}{2}\end{bmatrix}</div>

<h2>Substitution</h2><p>SymPy provides a means to substitute values in for the symbolic expressions. The specification requires an expression, a variable in the expression to substitute in for, and a new value. For example, this is one way to make a polynomial in a new variable:</p><pre class="sourceCode julia">@syms x y
ex = x^2 + 2x + 1
subs(ex, x, y)</pre>
<div class="well well-sm">
$$y^{2} + 2 y + 1$$</div>

<p>Substitution can also be numeric:</p><pre class="sourceCode julia">subs(ex, x, 0)</pre>
<div class="well well-sm">
$$1$$</div>

<p>The output has no free variables, but is still symbolic.</p><p>Expressions with more than one variable can have multiple substitutions, where each is expressed as a tuple:</p><pre class="sourceCode julia">x,y,z = symbols("x,y,z")
ex = x + y + z
subs(ex, (x,1), (y,pi))      </pre>
<div class="well well-sm">
$$z + 1 + \pi$$</div>

<p>Pairs can be used for substitution with:</p><pre class="sourceCode julia">subs(ex, x=>1, y=>pi)</pre>
<div class="well well-sm">
$$z + 1 + \pi$$</div>

<p>And, perhaps more conveniently, symbolic objects have their <code>call</code> method overloaded to allow substitution:</p><pre class="sourceCode julia">ex(x=>1, y=>pi)</pre>
<div class="well well-sm">
$$z + 1 + \pi$$</div>

<p>A straight call is also possble, where the order of the variables is determined by <code>free_symbols</code>:</p><pre class="sourceCode julia">ex(1, pi)</pre>
<div class="well well-sm">
$$y + 1 + \pi$$</div>

<p>When using the pipeline operator, <code>|&gt;</code>, is convenient, there is a curried form that allows the expression to be implicit:</p><pre class="sourceCode julia">ex |> subs(x, 1)</pre>
<div class="well well-sm">
$$y + z + 1$$</div>

<p>As <code>subs</code> is very similar in spirit to <code>Julia</code>'s <code>replace</code> function, that alias is provided:</p><pre class="sourceCode julia">ex |> replace(y, pi)</pre>
<div class="well well-sm">
$$x + z + 3.14159265358979$$</div>

<h2>Conversion from symbolic to numeric</h2><p>SymPy provides two identical means to convert a symbolic math expression to a number. One is <code>evalf</code>, the other <code>N</code>. Within <code>Julia</code> we decouple this, using <code>N</code> to also convert to a <code>Julian</code> value and <code>evalf</code> to leave the conversion as a symbolic object.  The <code>N</code> function converts symbolic integers, rationals, irrationals, and complex values, while attempting to find an appropriate <code>Julia</code> type for the value.</p><p>To see the difference, we use both on <code>PI</code>:</p><pre class="sourceCode julia">N(PI)  # floating-point value</pre>
<pre class="output">
3.141592653589793</pre>

<p>Whereas, while this may look the same, it is still symbolic:</p><pre class="sourceCode julia">evalf(PI)</pre>
<div class="well well-sm">
$$3.14159265358979$$</div>

<p>Both <code>N</code> and <code>evalf</code> allow for a precision argument to be passed through the second argument. This is how 30 digits of $\pi$ can be extracted:</p><pre class="sourceCode julia">N(PI, 30)</pre>
<pre class="output">
3.1415926535897932384626433832793</pre>

<p>Here <code>N</code> produces a <code>BigFloat</code> with a precision to match (basically) the specified number of digits. Whereas</p><pre class="sourceCode julia">evalf(PI, 30)</pre>
<div class="well well-sm">
$$3.14159265358979323846264338328$$</div>

<p>leaves the value as a symbolic object with 30 digits of accuracy.</p><p>Explicit conversion via <code>convert&#40;T, ex&#41;</code> can also be done, and is necessary at times if <code>N</code> does not give the desired type.</p><h2>Algebraic expressions</h2><p><code>SymPy</code> overloads many of <code>Julia</code>'s functions to work with symbolic objects, such as seen above with <code>asin</code>. The usual mathematical operations such as <code>&#43;</code>, <code>*</code>, <code>-</code>, <code>/</code> etc. work through <code>Julia</code>'s promotion mechanism, where numbers are promoted to symbolic objects, others dispatch internally to related <code>SymPy</code> functions.</p><p>In most all  cases, thinking about this distinction between numbers and symbolic numbers is unnecessary, as numeric values passed to <code>SymPy</code> functions are typically promoted to symbolic expressions. This conversion will take math constants to their corresponding <code>SymPy</code> counterpart, rational expressions to rational expressions, and floating point values to floating point values. However there are edge cases. An expression like <code>1//2 * pi * x</code> will differ from the seemingly identical  <code>1//2 * &#40;pi * x&#41;</code>. The former will produce a floating point value from <code>1//2 * pi</code> before being promoted to a symbolic instance. Using the symbolic value <code>PI</code> makes this expression work either way.</p><p>Most of <code>Julia</code>'s <a href="http://julia.readthedocs.org/en/latest/manual/mathematical-operations/#elementary-functions">mathematical</a> functions are overloaded to work with symbolic expressions. <code>Julia</code>'s generic definitions are used, as possible. This also introduces some edge cases. For example, <code>x^&#40;-2&#41;</code> will balk due to the negative, integer exponent, but either <code>x^&#40;-2//1&#41;</code> or <code>x^Sym&#40;-2&#41;</code> will work as expected, as the former call first dispatches to a generic defintion, but the latter two expressions do not.</p><p><code>SymPy</code> makes it very easy to work with polynomial and rational expressions. First we create some variables:</p><pre class="sourceCode julia">@syms x y z</pre>
<pre class="output">
(x,y,z)</pre>

<h3>The expand, factor, collect, and simplify functions</h3><p>A typical polynomial expression in a single variable can be written in two common ways, expanded or factored form. Using <code>factor</code> and <code>expand</code> can move between the two.</p><p>For example,</p><pre class="sourceCode julia">p = x^2 + 3x + 2
factor(p)</pre>
<div class="well well-sm">
$$\left(x + 1\right) \left(x + 2\right)$$</div>

<p>Or</p><pre class="sourceCode julia">expand(prod([(x-i) for i in 1:5]))</pre>
<div class="well well-sm">
$$x^{5} - 15 x^{4} + 85 x^{3} - 225 x^{2} + 274 x - 120$$</div>

<p>The <code>factor</code> function factors over the rational numbers, so something like this with obvious factors is not finished:</p><pre class="sourceCode julia">factor(x^2 - 2)</pre>
<div class="well well-sm">
$$x^{2} - 2$$</div>

<p>When expressions involve one or more variables, it can be convenient to be able to manipulate them. For example, if we define <code>q</code> by:</p><pre class="sourceCode julia">q = x*y + x*y^2 + x^2*y + x</pre>
<div class="well well-sm">
$$x^{2} y + x y^{2} + x y + x$$</div>

<p>Then we can collect the terms by the variable <code>x</code>:</p><pre class="sourceCode julia">collect(q, x)</pre>
<div class="well well-sm">
$$x^{2} y + x \left(y^{2} + y + 1\right)$$</div>

<p>or the variable <code>y</code>:</p><pre class="sourceCode julia">collect(q, y)</pre>
<div class="well well-sm">
$$x y^{2} + x + y \left(x^{2} + x\right)$$</div>

<p>These are identical expressions, though viewed differently. </p><p>A more broad-brush approach is to let <code>SymPy</code> simplify the values. In this case, the common value of <code>x</code> is factored out:</p><pre class="sourceCode julia">simplify(q)</pre>
<div class="well well-sm">
$$x \left(x y + y^{2} + y + 1\right)$$</div>

<p>The <code>simplify</code> function attempts to apply the dozens of functions related to simplification that are part of SymPy. It is also possible to apply these functions one at a time, for example <code>trigsimp</code> does trigonometric simplifications.</p><p>The SymPy tutorial illustrates that <code>expand</code> can also result in simplifications through this example:</p><pre class="sourceCode julia">expand((x + 1)*(x - 2) - (x - 1)*x)</pre>
<div class="well well-sm">
$$-2$$</div>

<p>These methods are not restricted to polynomial expressions and will work with other expressions. For example, <code>factor</code> identifies the following as a factorable object in terms of the variable <code>exp&#40;x&#41;</code>:</p><pre class="sourceCode julia">factor(exp(2x) + 3exp(x) + 2)</pre>
<div class="well well-sm">
$$\left(e^{x} + 1\right) \left(e^{x} + 2\right)$$</div>

<h2>Rational expressions: apart, together, cancel</h2><p>When working with rational expressions, SymPy does not do much simplification unless asked. For example this expression is not simplified:</p><pre class="sourceCode julia">r = 1/x + 1/x^2</pre>
<div class="well well-sm">
$$\frac{1}{x} + \frac{1}{x^{2}}$$</div>

<p>To put the terms of <code>r</code> over a common denominator, the <code>together</code> function is available:</p><pre class="sourceCode julia">together(r)</pre>
<div class="well well-sm">
$$\frac{1}{x^{2}} \left(x + 1\right)$$</div>

<p>The <code>apart</code> function does the reverse, creating a partial fraction decomposition from a ratio of polynomials:</p><pre class="sourceCode julia">apart( (4x^3 + 21x^2 + 10x + 12) /  (x^4 + 5x^3 + 5x^2 + 4x))</pre>
<div class="well well-sm">
$$\frac{2 x - 1}{x^{2} + x + 1} - \frac{1}{x + 4} + \frac{3}{x}$$</div>

<p>Some times SymPy will cancel factors, as here:</p><pre class="sourceCode julia">top = (x-1)*(x-2)*(x-3)
bottom = (x-1)*(x-4)
top/bottom</pre>
<div class="well well-sm">
$$\frac{\left(x - 3\right) \left(x - 2\right)}{x - 4}$$</div>

<p>(This might make math faculty a bit upset, but it is in line with student thinking.)</p><p>However, with expanded terms, the common factor of <code>&#40;x-1&#41;</code> is not cancelled:</p><pre class="sourceCode julia">r = expand(top) / expand(bottom)</pre>
<div class="well well-sm">
$$\frac{x^{3} - 6 x^{2} + 11 x - 6}{x^{2} - 5 x + 4}$$</div>

<p>The <code>cancel</code> function instructs SymPy to perform cancellations. It takes rational functions and puts them in a canonical $p/q$ form with no common (rational) factors and leading terms which are integers:</p><pre class="sourceCode julia">cancel(r)</pre>
<div class="well well-sm">
$$\frac{x^{2} - 5 x + 6}{x - 4}$$</div>

<h2>Powers</h2><p>The SymPy <a href="http://docs.sympy.org/dev/tutorial/simplification.html#powers">tutorial</a> offers a thorough explanation on powers and which get simplified and under what conditions. Basically</p><ul>
<li>$x^a x^b = x^{a+b}$
<p>is always true. However</p>
</li>
<li>$x^a y^a=(xy)^a$
<p>is only true with assumptions, such as $x,y \geq 0$ and $a$ is real, but not in general. For example, $x=y=-1$ and $a=1/2$ has $x^a \cdot y^a = i \cdot i =  -1$, where as $(xy)^a = 1$.</p>
</li>
<li>$(x^a)^b = x^{ab}$
<p>is only true with assumptions. For example $x=-1, a=2$, and $b=1/2$ gives $(x^a)^b = 1^{1/2} = 1$, whereas $x^{ab} = -1^1 = -1$.</p>
</li>
</ul><p>We see that with assumptions, the following expression does simplify to $0$:</p><pre class="sourceCode julia">@syms x y nonnegative=true a real=true
simplify(x^a * y^a - (x*y)^a)</pre>
<div class="well well-sm">
$$0$$</div>

<p>However, without assumptions this is not the case</p><pre class="sourceCode julia">x,y,a = symbols("x,y,a")
simplify(x^a * y^a - (x*y)^a)</pre>
<div class="well well-sm">
$$x^{a} y^{a} - \left(x y\right)^{a}$$</div>

<p>The <code>simplify</code> function calls <code>powsimp</code> to simplify powers, as above. The <code>powsimp</code> function has the keyword argument <code>force&#61;true</code> to force simplification even if assumptions are not specified:</p><pre class="sourceCode julia">powsimp(x^a * y^a - (x*y)^a, force=true)</pre>
<div class="well well-sm">
$$0$$</div>

<h2>Trigonometric simplification</h2><p>For trigonometric expressions, <code>simplify</code> will use <code>trigsimp</code> to simplify:</p><pre class="sourceCode julia">theta = symbols("theta", real=true)
p = cos(theta)^2 + sin(theta)^2</pre>
<div class="well well-sm">
$$\sin^{2}{\left (\theta \right )} + \cos^{2}{\left (\theta \right )}$$</div>

<p>Calling either <code>simplify</code> or <code>trigsimp</code> will apply the Pythagorean identity:</p><pre class="sourceCode julia">simplify(p)</pre>
<div class="well well-sm">
$$1$$</div>

<p>While often forgotten,  the <code>trigsimp</code> function is, of course,  aware of the double angle formulas:</p><pre class="sourceCode julia">simplify(sin(2theta) - 2sin(theta)*cos(theta))</pre>
<div class="well well-sm">
$$0$$</div>

<p>The <code>expand_trig</code> function will expand such expressions:</p><pre class="sourceCode julia">expand_trig(sin(2theta))</pre>
<div class="well well-sm">
$$2 \sin{\left (\theta \right )} \cos{\left (\theta \right )}$$</div>

<h2>Coefficients</h2><p>Returning to polynomials, there are a few functions to find various pieces of the polynomials. First we make a general quadratic polynomial:</p><pre class="sourceCode julia">a,b,c,x = symbols("a, b, c, x") 
p = a*x^2 + b*x + c</pre>
<div class="well well-sm">
$$a x^{2} + b x + c$$</div>

<p>The <code>coeff&#40;ex, monom&#41;</code> function will return the corresponding coefficient of the monomial:</p><pre class="sourceCode julia">coeff(p, x^2) # a 
coeff(p, x)   # b</pre>
<div class="well well-sm">
$$b$$</div>

<p>The constant can be found through substitution:</p><pre class="sourceCode julia">p(x=>0)</pre>
<div class="well well-sm">
$$c$$</div>

<p>Though one could use some trick like this to find all the coefficients:</p><pre class="sourceCode julia">Sym[[coeff(p, x^i) for i in N(degree(p)):-1:1]; p(x=>0)]</pre>
<div class="well well-sm">
\begin{bmatrix}a\\b\\c\end{bmatrix}</div>

<p>that is cumbersome, at best. SymPy has a function <code>coeffs</code>, but it is defined for polynomial types, so will fail on <code>p</code>:</p><pre class="sourceCode julia">coeffs(p) # fails</pre>
<pre class="output">
KeyError("coeffs")
</pre>

<p>Polynomials are a special class in SymPy and must be constructed. The <code>Poly</code> constructor can be used. As there is more than one free variable in <code>p</code>, we specify the variable <code>x</code> below:</p><pre class="sourceCode julia">q = Poly(p, x)
coeffs(q)</pre>
<pre class="output">
3-element Array{Any,1}:
 a
 b
 c</pre>

<h2>Polynomial roots: solve, real_roots, polyroots, nroots</h2><p>SymPy provides functions to find the roots of a polynomial. In general, a polynomial with real coefficients of degree $n$ will have $n$ roots when multiplicities and complex roots are accounted for. The number of real roots is consequently between $0$ and $n$.</p><p>For a <em>univariate</em> polynomial expression (a single variable), the real roots, when available, are returned by <code>real_roots</code>. For example,</p><pre class="sourceCode julia">real_roots(x^2 - 2)</pre>
<div class="well well-sm">
\begin{bmatrix}- \sqrt{2}\\\sqrt{2}\end{bmatrix}</div>

<p>Unlike <code>factor</code> – which only factors over rational factors – <code>real_roots</code> finds the two irrational roots here. It is well known (the <a href="http://en.wikipedia.org/wiki/Abel&#37;E2&#37;80&#37;93Ruffini_theorem">Abel-Ruffini theorem</a>) that for degree 5 polynomials, or higher, it is not always possible to express the roots in terms of radicals. However, when the roots are rational <code>SymPy</code> can have success:</p><pre class="sourceCode julia">p = (x-3)^2*(x-2)*(x-1)*x*(x+1)*(x^2 + x + 1)
real_roots(p)</pre>
<div class="well well-sm">
\begin{bmatrix}-1\\0\\1\\2\\3\\3\end{bmatrix}</div>

<p>In this example, the degree of <code>p</code> is 8, but only the 6 real roots returned, the double root of $3$ is accounted for. The two complex roots of <code>x^2 &#43; x&#43; 1</code> are not considered by this function. The complete set  of distinct roots can be found with <code>solve</code>:</p><pre class="sourceCode julia">solve(p)</pre>
<div class="well well-sm">
\begin{bmatrix}-1\\0\\1\\2\\3\\- \frac{1}{2} - \frac{\sqrt{3} i}{2}\\- \frac{1}{2} + \frac{\sqrt{3} i}{2}\end{bmatrix}</div>

<p>This finds the complex roots, but does not account for the double root. The <code>roots</code> function of SymPy does.</p><p>This particular function is not exported (as it conflicts with the <code>roots</code> function from the <code>Polynomials</code> package) but we can still access it using <code>p&#91;:roots&#93;&#40;&#41;</code> or its alias <code>polyroots</code>.</p><blockquote>
<p>Indexing with a symbol. When a symbolic expression is indexed by a symbol it returns a function which maps to a corresponding SymPy function. For example, <code>p&#91;:roots&#93;&#40;args...&#41;</code> will call <code>roots&#40;p, args...&#41;</code> within SymPy. For methods of SymPy objects, the same is true, so if <code>roots</code> were a class method, then the call would resolve to <code>p.roots&#40;args...&#41;</code>.</p>
</blockquote><p>The output of calling <code>polyroots</code> will be a dictionary whose keys are the roots and values the multiplicity.</p><pre class="sourceCode julia">polyroots(p)</pre>
<div class="well well-sm">
\begin{equation*}\begin{cases}1 & \text{=>} &1\\- \frac{1}{2} - \frac{\sqrt{3} i}{2} & \text{=>} &1\\3 & \text{=>} &2\\0 & \text{=>} &1\\-1 & \text{=>} &1\\- \frac{1}{2} + \frac{\sqrt{3} i}{2} & \text{=>} &1\\2 & \text{=>} &1\\\end{cases}\end{equation*}</div>

<p>When exact answers are not provided, the <code>polyroots</code> call is contentless:</p><pre class="sourceCode julia">p = x^5 - x + 1
polyroots(p)</pre>
<div class="well well-sm">
\begin{equation*}\begin{cases}\end{cases}\end{equation*}</div>

<p>Calling <code>solve</code> seems to produce very little as well:</p><pre class="sourceCode julia">rts = solve(p)</pre>
<div class="well well-sm">
\begin{bmatrix}\operatorname{CRootOf} {\left(x^{5} - x + 1, 0\right)}\\\operatorname{CRootOf} {\left(x^{5} - x + 1, 1\right)}\\\operatorname{CRootOf} {\left(x^{5} - x + 1, 2\right)}\\\operatorname{CRootOf} {\left(x^{5} - x + 1, 3\right)}\\\operatorname{CRootOf} {\left(x^{5} - x + 1, 4\right)}\end{bmatrix}</div>

<p>But in fact, <code>rts</code> contains lots of information. We can extract numeric values quite easily with <code>N</code>:</p><pre class="sourceCode julia">[N(r) for r in rts]     # or map(N, rts)</pre>
<pre class="output">
5-element Array{Number,1}:
         -1.1673     
 -0.181232-1.08395im 
 -0.181232+1.08395im 
  0.764884-0.352472im
  0.764884+0.352472im</pre>

<p>These are numeric approximations to irrational values. For numeric approximations to polynomial roots, the <code>nroots</code> function is also provided, though with this call the answers are still symbolic:</p><pre class="sourceCode julia">nroots(p)</pre>
<div class="well well-sm">
\begin{bmatrix}-1.16730397826142\\-0.181232444469875 - 1.08395410131771 i\\-0.181232444469875 + 1.08395410131771 i\\0.764884433600585 - 0.352471546031726 i\\0.764884433600585 + 0.352471546031726 i\end{bmatrix}</div>

<h2>The solve function</h2><p>The <code>solve</code> function is more general purpose than just finding roots of univariate polynomials. The function tries to solve for when an expression is 0, or a set of expressions are all 0.</p><p>For example, it can be used to solve when $\cos(x) = \sin(x)$:</p><pre class="sourceCode julia">solve(cos(x) - sin(x))</pre>
<div class="well well-sm">
\begin{bmatrix}- \frac{3 \pi}{4}\\\frac{\pi}{4}\end{bmatrix}</div>

<p>Though there are infinitely many correct solutions, these are within a certain range.</p><p>The <a href="http://docs.sympy.org/latest/modules/solvers/solveset.html">solveset</a> function appears in version 1.0 of SymPy and is an intended replacement for <code>solve</code>. Here we see it gives all solutions:</p><pre class="sourceCode julia">u = solveset(cos(x) - sin(x))</pre>
<div class="well well-sm">
$$\left\{2 n \pi - \frac{3 \pi}{4}\; |\; n \in \mathbb{Z}\right\} \cup \left\{2 n \pi + \frac{\pi}{4}\; |\; n \in \mathbb{Z}\right\}$$</div>

<p>The output of <code>solveset</code> is a set, rather than a vector or dictionary. To get the values requires some work. For <em>finite sets</em> we collect the elements with <code>elements</code>:</p><pre class="sourceCode julia">v = solveset(x^2 - 4)
elements(v)       </pre>
<pre class="output">
2-element Array{Any,1}:
 -2
  2</pre>

<p>The <code>elements</code> function does not work for more complicated sets, such as <code>u</code>. For these, the <code>contains</code> method may be useful.</p><p>Solving within Sympy has limits. For example, there is no symbolic solution here:</p><pre class="sourceCode julia">solve(cos(x) - x)</pre>
<pre class="output">
PyError (:PyObject_Call) <type 'exceptions.NotImplementedError'>
NotImplementedError('multiple generators [x, cos(x)]\nNo algorithms are implemented to solve equation -x + cos(x)',)
  File "/Users/verzani/.julia/v0.5/Conda/deps/usr/lib/python2.7/site-packages/sympy/solvers/solvers.py", line 1053, in solve
    solution = _solve(f[0], *symbols, **flags)
  File "/Users/verzani/.julia/v0.5/Conda/deps/usr/lib/python2.7/site-packages/sympy/solvers/solvers.py", line 1619, in _solve
    raise NotImplementedError('\n'.join([msg, not_impl_msg % f]))

</pre>

<p>For such, a numeric method would be needed, say:</p><pre class="sourceCode julia">nsolve(cos(x) - x, 1)</pre>
<pre class="output">
7.390851332151606416553120876738734040134117589007574649656806357732846548834992e-01</pre>

<p>Though it can't solve everything, the <code>solve</code> function can also solve equations of a more general type. For example, here it is used to derive the quadratic equation:</p><pre class="sourceCode julia">a,b,c  = symbols("a,b,c", real=true)
p = a*x^2 + b*x + c
solve(p, x)</pre>
<div class="well well-sm">
\begin{bmatrix}\frac{1}{2 a} \left(- b + \sqrt{- 4 a c + b^{2}}\right)\\- \frac{1}{2 a} \left(b + \sqrt{- 4 a c + b^{2}}\right)\end{bmatrix}</div>

<p>The extra argument <code>x</code> is passed to <code>solve</code> so that <code>solve</code> knows which variable to solve for.</p><p>The <code>solveset</code> function is similar:</p><pre class="sourceCode julia">solveset(p, x)</pre>
<div class="well well-sm">
$$\left\{- \frac{b}{2 a} - \frac{1}{2 a} \sqrt{- 4 a c + b^{2}}, - \frac{b}{2 a} + \frac{1}{2 a} \sqrt{- 4 a c + b^{2}}\right\}$$</div>

<p>If the <code>x</code> value is not given, <code>solveset</code> will complain and  <code>solve</code> tries to find a solution with all the free variables:</p><pre class="sourceCode julia">solve(p)</pre>
<pre class="output">
1-element Array{Dict{SymPy.Sym,SymPy.Sym},1}:
 Dict(a=>-(b*x + c)/x^2)</pre>

<p>Systems of equations can be solved as well. We specify them within a vector of expressions, <code>&#91;ex1, ex2, ..., exn&#93;</code> where a found solution is one where all the expressions are 0. For example, to solve this linear system: $2x + 3y = 6, 3x - 4y=12$, we have:</p><pre class="sourceCode julia">x, y = symbols("x,y", real=true)
exs = [2x+3y-6, 3x-4y-12]
d = solve(exs)</pre>
<div class="well well-sm">
\begin{equation*}\begin{cases}y & \text{=>} &- \frac{6}{17}\\x & \text{=>} &\frac{60}{17}\\\end{cases}\end{equation*}</div>

<p>We can "check our work" by plugging into each equation. We take advantage of how the <code>subs</code> function allows us to pass in a dictionary:</p><pre class="sourceCode julia">map(ex -> subs(ex, d), exs)</pre>
<div class="well well-sm">
\begin{bmatrix}0\\0\end{bmatrix}</div>

<p>In the previous example, the system had two equations and two unknowns. When that is not the case, one can specify the variables to solve for as a vector. In this example, we find a quadratic polynomial that approximates $\cos(x)$ near $0$:</p><pre class="sourceCode julia">a,b,c,h = symbols("a,b,c,h", real=true)
p = a*x^2 + b*x + c
fn = cos
exs = [fn(0*h)-p(x=>0), fn(h)-p(x => h), fn(2h)-p(x => 2h)]
d = solve(exs, [a,b,c])</pre>
<div class="well well-sm">
\begin{equation*}\begin{cases}a & \text{=>} &\frac{1}{2 h^{2}} \left(- 2 \cos{\left (h \right )} + \cos{\left (2 h \right )} + 1\right)\\b & \text{=>} &\frac{1}{2 h} \left(4 \cos{\left (h \right )} - \cos{\left (2 h \right )} - 3\right)\\c & \text{=>} &1\\\end{cases}\end{equation*}</div>

<p>Again, a dictionary is returned. The polynomial itself can be found by substituting back in for <code>a</code>, <code>b</code>, and <code>c</code>:</p><pre class="sourceCode julia">quad_approx = subs(p, d)</pre>
<div class="well well-sm">
$$1 + \frac{x}{2 h} \left(4 \cos{\left (h \right )} - \cos{\left (2 h \right )} - 3\right) + \frac{x^{2}}{2 h^{2}} \left(- 2 \cos{\left (h \right )} + \cos{\left (2 h \right )} + 1\right)$$</div>

<p>(Taking the limit as $h$ goes to 0 produces the answer $1 - x^2/2$.)</p><p>Finally for <code>solve</code>, we show one way to re-express the polynomial $a_2x^2 + a_1x + a_0$ as $b_2(x-c)^2 + b_1(x-c) + b_0$ using <code>solve</code> (and not, say, an expansion theorem.)</p><pre class="sourceCode julia">n = 3
x, c = symbols("x,c")
as = Sym["a$i" for i in 0:(n-1)]
bs = Sym["b$i" for i in 0:(n-1)]
p = sum([as[i+1]*x^i for i in 0:(n-1)])
q = sum([bs[i+1]*(x-c)^i for i in 0:(n-1)])
solve(p-q, bs)</pre>
<div class="well well-sm">
\begin{equation*}\begin{cases}b_{0} & \text{=>} &a_{0} + a_{1} c + a_{2} c^{2}\\b_{1} & \text{=>} &a_{1} + 2 a_{2} c\\b_{2} & \text{=>} &a_{2}\\\end{cases}\end{equation*}</div>

<h3>Solving using logical operators</h3><p>The <code>solve</code> function does not need to just solve <code>ex &#61; 0</code>. There are other means to specify an equation. Ideally, it would be nice to say <code>ex1 &#61;&#61; ex2</code>, but the interpretation of <code>&#61;&#61;</code> is not for this. Rather, <code>SymPy</code> introduces <code>Eq</code> for equality. So this expression</p><pre class="sourceCode julia">solve(Eq(x, 1))</pre>
<div class="well well-sm">
\begin{bmatrix}1\end{bmatrix}</div>

<p>gives 1, as expected from solving <code>x &#61;&#61; 1</code>.</p><p>In addition to <code>Eq</code>, there are <code>Lt</code>, <code>Le</code>, <code>Ge</code>, <code>Gt</code>. The Unicode operators are not aliased to these, but there are alternatives <code>\ll&#91;tab&#93;</code>, <code>\leqq&#91;tab&#93;</code>, <code>\Equal&#91;tab&#93;</code>, <code>\geqq&#91;tab&#93;</code>, <code>\gg&#91;tab&#93;</code> and <code>\neg&#91;tab&#93;</code> to negate.</p><p>So, the above could have been written with the following nearly identical expression, though it is entered with <code>\Equal&#91;tab&#93;</code>. </p><pre class="sourceCode julia">solve(x ⩵ 1)</pre>
<div class="well well-sm">
\begin{bmatrix}1\end{bmatrix}</div>

<p>Here is an alternative way of asking a previous question on a pair of linear equations:</p><pre class="sourceCode julia">x, y = symbols("x,y", real=true)
exs = [2x+3y ⩵ 6, 3x-4y ⩵ 12]    ## Using \Equal[tab]
d = solve(exs)</pre>
<div class="well well-sm">
\begin{equation*}\begin{cases}y & \text{=>} &- \frac{6}{17}\\x & \text{=>} &\frac{60}{17}\\\end{cases}\end{equation*}</div>

<h2>Plotting</h2><p>The <code>Plots</code> package allows many 2-dimensional plots of <code>SymPy</code> objects to be agnostic as to a backend plotting package.  <code>SymPy</code> provides recipes that allow symbolic expressions to be used where functions are part of the <code>Plots</code> interface. [See the help page for <code>sympy_plotting</code>.]</p><p>In particular, the following methods of <code>plot</code> are defined:</p><ul>
<li><p><code>plot&#40;ex::Sym, a, b&#41;</code> will plot the expression of single variable over the interval <code>&#91;a,b&#93;</code></p>
</li>
<li><p><code>plot&#33;&#40;ex::Sym, a, b&#41;</code> will add to the current plot a plot of  the expression of single variable over the interval <code>&#91;a,b&#93;</code></p>
</li>
<li><p><code>plot&#40;exs::Vector&#123;Sym&#125;, a, b&#41;</code> will plot each expression over <code>&#91;a,b&#93;</code></p>
</li>
<li><p><code>plot&#40;ex1, ex2, a, b&#41;</code> will plot a parametric plot of the two expressions over the interval <code>&#91;a,b&#93;</code>.</p>
</li>
<li><p><code>contour&#40;xs, ys, ex::Sym&#41;</code> will make a contour plot of the expression of two variables over the grid specifed by the <code>xs</code> and <code>ys</code>.</p>
</li>
<li><p><code>surface&#40;xs, ys, ex::Sym&#41;</code> will make a surface plot of the expression of two variables over the grid specifed by the <code>xs</code> and <code>ys</code>.</p>
</li>
</ul><p>For example:</p><pre class="sourceCode julia">x = symbols("x")
using Plots
#
plot(x^2 - 2, -2,2)</pre>
    <div id="efb51c59-1b68-47d6-bcdc-0fbde3fe26a2" style="width:600px;height:400px;"></div>
    <script>
    PLOT = document.getElementById('efb51c59-1b68-47d6-bcdc-0fbde3fe26a2');
    Plotly.plot(PLOT, [{"yaxis":"y","y":[2.0,1.8400163248648096,1.683297622691562,1.5298438934802574,1.3796551372308952,1.2327313539434752,1.089072543617998,0.9486787062544635,0.8115498418528717,0.6776859504132235,0.5470870319355168,0.4197530864197532,0.29568411386593185,0.1748801142740537,0.05734108764411783,-0.056932966023874876,-0.1679420467299253,-0.2756861544740332,-0.3801652892561984,-0.4813794510764209,-0.5793286399347009,-0.6740128558310374,-0.765432098765432,-0.853586368737884,-0.938475665748393,-1.0200999897969594,-1.0984593408835832,-1.1735537190082646,-1.2453831241710027,-1.3139475563717988,-1.379247015610652,-1.4412815018875624,-1.5000510152025304,-1.5555555555555556,-1.607795122946638,-1.656769717375778,-1.7024793388429753,-1.7449239873482298,-1.7841036628915417,-1.820018365472911,-1.8526680950923375,-1.8820528517498214,-1.9081726354453628,-1.9310274461789614,-1.9506172839506173,-1.9669421487603307,-1.9800020406081011,-1.989796959493929,-1.9963269054178145,-1.9995918783797573,-1.9995918783797573,-1.9963269054178145,-1.989796959493929,-1.9800020406081011,-1.9669421487603307,-1.9506172839506173,-1.9310274461789614,-1.9081726354453628,-1.8820528517498214,-1.8526680950923375,-1.820018365472911,-1.7841036628915417,-1.7449239873482298,-1.7024793388429753,-1.656769717375778,-1.607795122946638,-1.5555555555555556,-1.5000510152025304,-1.4412815018875624,-1.379247015610652,-1.3139475563717988,-1.2453831241710027,-1.1735537190082646,-1.0984593408835832,-1.0200999897969594,-0.938475665748393,-0.853586368737884,-0.765432098765432,-0.6740128558310374,-0.5793286399347009,-0.4813794510764209,-0.3801652892561984,-0.2756861544740332,-0.1679420467299253,-0.056932966023874876,0.05734108764411783,0.1748801142740537,0.29568411386593185,0.4197530864197532,0.5470870319355168,0.6776859504132235,0.8115498418528717,0.9486787062544635,1.089072543617998,1.2327313539434752,1.3796551372308952,1.5298438934802574,1.683297622691562,1.8400163248648096,2.0],"showlegend":true,"name":"y1","type":"scatter","xaxis":"x","line":{"width":1,"dash":"solid","color":"rgba(0, 154, 250, 1.000)","shape":"linear"},"x":[-2.0,-1.9595959595959596,-1.9191919191919191,-1.878787878787879,-1.8383838383838385,-1.797979797979798,-1.7575757575757576,-1.7171717171717171,-1.6767676767676767,-1.6363636363636365,-1.595959595959596,-1.5555555555555556,-1.5151515151515151,-1.4747474747474747,-1.4343434343434343,-1.393939393939394,-1.3535353535353536,-1.3131313131313131,-1.2727272727272727,-1.2323232323232323,-1.1919191919191918,-1.1515151515151516,-1.1111111111111112,-1.0707070707070707,-1.0303030303030303,-0.98989898989899,-0.9494949494949495,-0.9090909090909091,-0.8686868686868687,-0.8282828282828283,-0.7878787878787878,-0.7474747474747475,-0.7070707070707071,-0.6666666666666666,-0.6262626262626263,-0.5858585858585859,-0.5454545454545454,-0.5050505050505051,-0.46464646464646464,-0.42424242424242425,-0.3838383838383838,-0.3434343434343434,-0.30303030303030304,-0.26262626262626265,-0.2222222222222222,-0.18181818181818182,-0.1414141414141414,-0.10101010101010101,-0.06060606060606061,-0.020202020202020204,0.020202020202020204,0.06060606060606061,0.10101010101010101,0.1414141414141414,0.18181818181818182,0.2222222222222222,0.26262626262626265,0.30303030303030304,0.3434343434343434,0.3838383838383838,0.42424242424242425,0.46464646464646464,0.5050505050505051,0.5454545454545454,0.5858585858585859,0.6262626262626263,0.6666666666666666,0.7070707070707071,0.7474747474747475,0.7878787878787878,0.8282828282828283,0.8686868686868687,0.9090909090909091,0.9494949494949495,0.98989898989899,1.0303030303030303,1.0707070707070707,1.1111111111111112,1.1515151515151516,1.1919191919191918,1.2323232323232323,1.2727272727272727,1.3131313131313131,1.3535353535353536,1.393939393939394,1.4343434343434343,1.4747474747474747,1.5151515151515151,1.5555555555555556,1.595959595959596,1.6363636363636365,1.6767676767676767,1.7171717171717171,1.7575757575757576,1.797979797979798,1.8383838383838385,1.878787878787879,1.9191919191919191,1.9595959595959596,2.0],"mode":"lines"}], {"yaxis":{"type":"-","titlefont":{"size":15,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"title":"","tickfont":{"size":11,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"ticks":"inside","tickmode":"array","showgrid":true,"tickvals":[-1.0,0.0,1.0,2.0],"domain":[0.03762029746281716,0.9901574803149606],"ticktext":["-1","0","1","2"],"tickangle":0,"zeroline":false,"linecolor":"rgba(0, 0, 0, 1.000)","tickcolor":"rgba(0, 0, 0, 1.000)","anchor":"x"},"annotations":[],"width":600,"plot_bgcolor":"rgba(255, 255, 255, 1.000)","showlegend":true,"legend":{"bgcolor":"rgba(255, 255, 255, 1.000)","y":1.0,"font":{"size":11,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"bordercolor":"rgba(0, 0, 0, 1.000)","x":1.0},"xaxis":{"type":"-","titlefont":{"size":15,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"title":"","tickfont":{"size":11,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"ticks":"inside","tickmode":"array","showgrid":true,"tickvals":[-2.0,-1.0,0.0,1.0,2.0],"domain":[0.03619130941965587,0.9934383202099738],"ticktext":["-2","-1","0","1","2"],"tickangle":0,"zeroline":false,"linecolor":"rgba(0, 0, 0, 1.000)","tickcolor":"rgba(0, 0, 0, 1.000)","anchor":"y"},"paper_bgcolor":"rgba(255, 255, 255, 1.000)","height":400,"margin":{"r":0,"l":0,"b":0,"t":20}});
    </script>
<p>Or a parametric plot:</p><pre class="sourceCode julia">plot(sin(2x), cos(3x), 0, 4pi)</pre>
    <div id="14ad6034-269e-4e7c-b39b-e2c097b8d155" style="width:600px;height:400px;"></div>
    <script>
    PLOT = document.getElementById('14ad6034-269e-4e7c-b39b-e2c097b8d155');
    Plotly.plot(PLOT, [{"yaxis":"y","y":[1,0.9821093199149804,0.9290774325277306,0.8428018909013506,0.7263697512646394,0.5839471139413063,0.42063005461378417,0.24226227980378318,0.05522603110450816,-0.13378628010447902,-0.3180115362392381,-0.4908579070575935,-0.6461407143112099,-0.77830372794553,-0.8826179755685469,-0.9553509515151947,-0.9939001709768878,-0.9968862904477932,-0.9642024625116118,-0.8970181589874636,-0.7977373256375195,-0.6699123657178551,-0.5181170301580778,-0.3477827625319824,-0.1650043546187989,0.023678133536623857,0.21151338586781926,0.3917804015584929,0.5580289815934403,0.7043105256526725,0.8253908811219757,0.9169376281927886,0.9756750997357736,0.9995015891261737,0.9875645521655236,0.9402911122726758,0.8593727774269122,0.747704915711709,0.6092831551065171,0.449060414482921,0.27276968143060365,0.08671887816355096,-0.10243484451661349,-0.2879233091311653,-0.4631094861203471,-0.6217249757884948,-0.7580943001712455,-0.8673379793567147,-0.9455471259136673,-0.989923310200557,-0.9988786919844436,-0.972092835524257,-0.910524175197462,-0.8163757214143996,-0.6930162339093321,-0.5448596829350706,-0.3772073114035758,-0.19605794920397915,-0.007893366909712247,0.18055365078890184,0.36254021327862285,0.5315545938209018,0.6815492279916331,0.807157103761985,0.9038838004888234,0.9682683053985073,0.9980068533314933,0.9920353585932581,0.9505674893877825,0.8750870224785937,0.7682947516379717,0.6340118495722394,0.47704314116489016,0.3030051803068747,0.11812528195890788,-0.07098129964801639,-0.25754807380689576,-0.4348994275757909,-0.5966894882888544,-0.7371291875117791,-0.851193401784494,-0.934800758373598,-0.9849596723401104,-0.9998753895176573,-0.9790142052577144,-0.9231225611078604,-0.8342003361179189,-0.7154292884473727,-0.5710592077306961,-0.40625585182379015,-0.22691610896159098,-0.03945699907625321,0.14941393590426086,0.3329386370297614,0.504550340869178,0.6581085472380385,0.7881187346471933,0.8899289617551795,0.9598963201568567,0.9955172826011062,0.995517282601106,0.959896320156857,0.8899289617551802,0.7881187346471942,0.6581085472380369,0.5045503408691763,0.3329386370297628,0.14941393590426233,-0.03945699907625174,-0.22691610896158954,-0.40625585182378876,-0.5710592077306949,-0.7154292884473716,-0.834200336117918,-0.9231225611078598,-0.9790142052577149,-0.9998753895176573,-0.9849596723401106,-0.9348007583735986,-0.8511934017844948,-0.7371291875117814,-0.5966894882888556,-0.43489942757579225,-0.2575480738068989,-0.07098129964801786,0.11812528195890641,0.3030051803068716,0.4770431411648889,0.6340118495722383,0.7682947516379685,0.8750870224785938,0.9505674893877831,0.9920353585932575,0.9980068533314933,0.9682683053985085,0.9038838004888241,0.8071571037619838,0.6815492279916368,0.5315545938209015,0.3625402132786275,0.18055365078890503,-0.007893366909710778,-0.19605794920397596,-0.3772073114035744,-0.5448596829350694,-0.6930162339093311,-0.8163757214143967,-0.9105241751974621,-0.9720928355242571,-0.9988786919844435,-0.989923310200557,-0.945547125913669,-0.8673379793567155,-0.7580943001712441,-0.6217249757884988,-0.4631094861203484,-0.28792330913117015,-0.10243484451661672,0.08671887816355127,0.27276968143060054,0.44906041448291806,0.6092831551065152,0.747704915711708,0.8593727774269095,0.9402911122726753,0.9875645521655225,0.9995015891261738,0.9756750997357735,0.9169376281927885,0.8253908811219756,0.704310525652671,0.5580289815934445,0.3917804015584975,0.2115133858678233,0.02367813353662799,-0.1650043546187957,-0.3477827625319794,-0.518117030158075,-0.6699123657178534,-0.7977373256375181,-0.8970181589874627,-0.9642024625116113,-0.9968862904477931,-0.9939001709768887,-0.9553509515151949,-0.8826179755685502,-0.7783037279455296,-0.646140714311209,-0.49085790705759247,-0.3180115362392433,-0.13378628010447696,0.055226031104503595,0.24226227980378603,0.42063005461378083,0.5839471139413035,0.7263697512646373,0.8428018909013493,0.9290774325277297,0.9821093199149814,1.0],"showlegend":true,"name":"y1","type":"scatter","xaxis":"x","line":{"width":1,"dash":"solid","color":"rgba(0, 154, 250, 1.000)","shape":"linear"},"x":[0,0.12595970506771756,0.24991296237030836,0.3698852854165468,0.48396560114283876,0.5903366919365284,0.6873041426091843,0.773323331215339,0.8470240337722988,0.9072322509454813,0.9529889090158392,0.9835651372363698,0.9984738782203788,0.9974776464163387,0.9805923110824041,0.9480868435005095,0.9004790324567516,0.8385272362373773,0.7632183025251686,0.6757518486236088,0.5775211524135886,0.47009095843600296,0.3551725526334286,0.23459650279236888,0.11028349891127508,-0.01578624201363702,-0.14160451942495178,-0.2651671383986785,-0.38450583503201075,-0.49771962947568316,-0.6030051071796143,-0.6986851459933063,-0.7832356315188904,-0.8553097351604116,-0.9137593681374367,-0.9576534697159295,-0.9862928383380026,-0.9992212694012756,-0.9962328222710067,-0.9773751007667073,-0.9429484948674371,-0.8935013957148743,-0.8298214601357258,-0.7529230638333771,-0.6640311431104317,-0.5645616825119181,-0.4560991592070162,-0.34037130340411265,-0.21922157684769145,-0.09457980779484539,0.03156854976481064,0.15721404296725056,0.2803552292170147,0.399030557732341,0.5113496156423264,0.6155232409081793,0.7098920223913328,0.7929527330827784,0.8633822754312234,0.9200587573381744,0.9620793630944627,0.9887747345870026,0.9997196336934779,0.9947397160206569,0.9739143080855378,0.9375751437008252,0.886301079693209,0.8209088751292632,0.7424401809292829,0.6521449471151867,0.5514615120031078,0.44199369050557974,0.32548522651021267,0.20379201629015345,0.07885254540747641,-0.04734298997155772,-0.17278438547409894,-0.29547344963466915,-0.4134558336521329,-0.5248521627644684,-0.6278879733408582,-0.720921978814716,-0.8024722142015774,-0.8712396427384599,-0.9261288486078413,-0.9662654861260218,-0.991010207442792,-0.9999688468941563,-0.9929986997786697,-0.9702107955409862,-0.9319681291524349,-0.8788798788614605,-0.8117917024210212,-0.7317722663670775,-0.6400962229271068,-0.5382239057242885,-0.42777806771020976,-0.3105180318741694,-0.18831166648971873,-0.06310563131267308,0.06310563131267387,0.18831166648971776,0.31051803187416843,0.4277780677102089,0.5382239057242891,0.6400962229271074,0.7317722663670768,0.8117917024210206,0.87887987886146,0.9319681291524344,0.9702107955409863,0.9929986997786696,0.9999688468941563,0.9910102074427922,0.966265486126022,0.926128848607841,0.8712396427384594,0.8024722142015781,0.7209219788147168,0.627887973340859,0.5248521627644693,0.4134558336521338,0.2954734496346701,0.1727843854740999,0.0473429899715587,-0.07885254540747542,-0.2037920162901525,-0.32548522651021006,-0.44199369050557885,-0.5514615120031054,-0.6521449471151859,-0.7424401809292834,-0.8209088751292617,-0.8863010796932086,-0.9375751437008243,-0.9739143080855376,-0.994739716020657,-0.999719633693478,-0.9887747345870026,-0.9620793630944634,-0.9200587573381748,-0.8633822754312248,-0.7929527330827795,-0.7098920223913329,-0.6155232409081808,-0.5113496156423272,-0.39903055773234347,-0.28035522921701567,-0.15721404296725064,-0.031568549764812506,0.0945798077948444,0.21922157684768875,0.3403713034041117,0.45609915920701616,0.5645616825119165,0.6640311431104309,0.7529230638333754,0.8298214601357252,0.8935013957148742,0.9429484948674365,0.9773751007667072,0.9962328222710065,0.9992212694012756,0.9862928383380032,0.95765346971593,0.9137593681374383,0.8553097351604113,0.783235631518891,0.6986851459933061,0.6030051071796159,0.4977196294756832,0.3845058350320129,0.26516713839867945,0.1416045194249545,0.015786242013638443,-0.11028349891127145,-0.23459650279236663,-0.35517255263342773,-0.4700909584360037,-0.5775211524135875,-0.6757518486236088,-0.7632183025251672,-0.838527236237377,-0.9004790324567504,-0.9480868435005091,-0.9805923110824034,-0.9974776464163386,-0.9984738782203788,-0.9835651372363697,-0.9529889090158397,-0.9072322509454812,-0.8470240337722998,-0.7733233312153394,-0.6873041426091862,-0.5903366919365293,-0.48396560114284154,-0.36988528541654847,-0.2499129623703121,-0.12595970506771642,-9.797174393178826e-16],"mode":"lines"}], {"yaxis":{"type":"-","titlefont":{"size":15,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"title":"","tickfont":{"size":11,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"ticks":"inside","tickmode":"array","showgrid":true,"tickvals":[-0.5,0.0,0.5,1.0],"domain":[0.03762029746281716,0.9901574803149606],"ticktext":["-0.5","0.0","0.5","1.0"],"tickangle":0,"zeroline":false,"linecolor":"rgba(0, 0, 0, 1.000)","tickcolor":"rgba(0, 0, 0, 1.000)","anchor":"x"},"annotations":[],"width":600,"plot_bgcolor":"rgba(255, 255, 255, 1.000)","showlegend":true,"legend":{"bgcolor":"rgba(255, 255, 255, 1.000)","y":1.0,"font":{"size":11,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"bordercolor":"rgba(0, 0, 0, 1.000)","x":1.0},"xaxis":{"type":"-","titlefont":{"size":15,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"title":"","tickfont":{"size":11,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"ticks":"inside","tickmode":"array","showgrid":true,"tickvals":[-0.5,0.0,0.5],"domain":[0.0658209390492855,0.9934383202099738],"ticktext":["-0.5","0.0","0.5"],"tickangle":0,"zeroline":false,"linecolor":"rgba(0, 0, 0, 1.000)","tickcolor":"rgba(0, 0, 0, 1.000)","anchor":"y"},"paper_bgcolor":"rgba(255, 255, 255, 1.000)","height":400,"margin":{"r":0,"l":0,"b":0,"t":20}});
    </script>
<p>For plotting with other plotting packages, it is generally faster to first call <code>lambdify</code> on the expression and then generate <code>y</code> values with the resulting <code>Julia</code> function.</p><hr /><p>In addition, with <code>PyPlot</code> a few other plotting functions from <code>SymPy</code> are available from its interface to <code>MatplotLib</code>:</p><ul>
<li><p><code>plot_parametric_surface&#40;ex1::Sym, ex2::Sym, ex3::Sym&#41;, &#40;uvar, a0, b0&#41;, &#40;vvar, a1, b1&#41;&#41;</code> – make a surface plot of the expressions parameterized by the region <code>&#91;a0,b0&#93; x &#91;a1,b1&#93;</code>. The default region is <code>&#91;-5,5&#93;x&#91;-5,5&#93;</code> where the ordering of the variables is given by <code>free_symbols&#40;ex&#41;</code>.</p>
</li>
<li><p><code>plot_implicit&#40;predictate, &#40;xvar, a0, b0&#41;, &#40;yvar, a1, b1&#41;&#41;</code> – make</p>
</li>
</ul><p>an implicit equation plot of the expressions over the region <code>&#91;a0,b0&#93; x &#91;a1,b1&#93;</code>. The default region is <code>&#91;-5,5&#93;x&#91;-5,5&#93;</code> where the ordering of the variables is given by <code>free_symbols&#40;ex&#41;</code>.  To create predicates from the variable, the functions <code>Lt</code>, <code>Le</code>, <code>Eq</code>, <code>Ge</code>, and <code>Gt</code> can be used, as with <code>Lt&#40;x*y, 1&#41;</code>. For infix notation, unicode operators can be used: <code>\ll&lt;tab&gt;</code>, <code>\leqq&lt;tab&gt;</code>, <code>\Equal&lt;tab&gt;</code>, <code>\geqq&lt;tab&gt;</code>, and <code>\gg&lt;tab&gt;</code>. For example, <code>x*y ≪ 1</code>.  To combine terms, the unicode <code>\vee&lt;tab&gt;</code> (for "or"), <code>\wedge&lt;tab&gt;</code> (for "and") can be used.</p><h2>Calculus</h2><p><code>SymPy</code> has many of the basic operations of calculus provided through a relatively small handful of functions.</p><h3>Limits</h3><p>Limits are computed by the <code>limit</code> function which takes an expression, a variable and a value, and optionally a direction specified by either <code>dir&#61;&quot;&#43;&quot;</code> or <code>dir&#61;&quot;-&quot;</code>.</p><p>For example, this shows Gauss was right:</p><pre class="sourceCode julia">limit(sin(x)/x, x, 0)</pre>
<div class="well well-sm">
$$1$$</div>

<p>Alternatively, the second and third arguments can be specified as a pair:</p><pre class="sourceCode julia">limit(sin(x)/x, x=>0)</pre>
<div class="well well-sm">
$$1$$</div>

<p>Limits at infinity are done by using <code>oo</code> for $\infty$:</p><pre class="sourceCode julia">limit((1+1/x)^x, x => oo)</pre>
<div class="well well-sm">
$$e$$</div>

<p>This example computes what L'Hopital reportedly paid a Bernoulli for</p><pre class="sourceCode julia">a = symbols("a", positive=true)
ex = (sqrt(2a^3*x-x^4) - a*(a^2*x)^(1//3)) / (a - (a*x^3)^(1//4))</pre>
<div class="well well-sm">
$$\frac{- a^{\frac{5}{3}} \root{3}{x} + \sqrt{2 a^{3} x - x^{4}}}{- \root{4}{a} \root{4}{x^{3}} + a}$$</div>

<p>Substituting $x=a$ gives an indeterminate form:</p><pre class="sourceCode julia">ex(x=>a)         # or subs(ex, x, a)</pre>
<div class="well well-sm">
$$\mathrm{NaN}$$</div>

<p>We can see it is of the form $0/0$:</p><pre class="sourceCode julia">subs(denom(ex), x, a), subs(numer(ex), x, a)</pre>
<pre class="output">
(0,0)</pre>

<p>And we get</p><pre class="sourceCode julia">limit(ex, x => a)</pre>
<div class="well well-sm">
$$\frac{16 a}{9}$$</div>

<p>In a previous example, we defined <code>quad_approx</code>:</p><pre class="sourceCode julia">quad_approx</pre>
<div class="well well-sm">
$$1 + \frac{x}{2 h} \left(4 \cos{\left (h \right )} - \cos{\left (2 h \right )} - 3\right) + \frac{x^{2}}{2 h^{2}} \left(- 2 \cos{\left (h \right )} + \cos{\left (2 h \right )} + 1\right)$$</div>

<p>The limit as <code>h</code> goes to $0$ gives <code>1 - x^2/2</code>, as expected:</p><pre class="sourceCode julia">limit(quad_approx, h => 0)</pre>
<div class="well well-sm">
$$- \frac{x^{2}}{2} + 1$$</div>

<h4>Left and right limits</h4><p>The limit is defined when both the left and right limits exist and are equal. But left and right limits can exist and not be equal. The <code>sign</code> function is $1$ for positive $x$, $-1$ for negative $x$ and $0$ when $x$ is 0. It should not have a limit at $0$:</p><pre class="sourceCode julia">limit(sign(x), x => 0)</pre>
<div class="well well-sm">
$$1$$</div>

<p>Oops. Well, the left and right limits are different anyways:</p><pre class="sourceCode julia">limit(sign(x), x => 0, dir="-"), limit(sign(x), x => 0, dir="+")</pre>
<pre class="output">
(-1,1)</pre>

<p>(The <code>limit</code> function finds the <em>right</em> limit by default. To be careful, either plot or check that both the left and right limit exist and are equal.)</p><h4>Operator interface</h4><p>For univariate functions there is an "operator" interface, where we pass a function object as the first argument and the value for <code>c</code> as the second (the variable is implicit, as <code>f</code> has only one).</p><pre class="sourceCode julia">f(x) = sin(5x)/x
limit(f, 0)</pre>
<pre class="output">
MethodError(*,(5,Symbol("##282")))
</pre>

<h4>Numeric limits</h4><p>The <code>limit</code> function uses the <a href="http://docs.sympy.org/latest/modules/series.html#the-gruntz-algorithm">Gruntz</a> algorithm. It is far more reliable then simple numeric attempts at limits. An example of Gruntz is the right limit at $0$ of the function:</p><pre class="sourceCode julia">f(x) = 1/x^(log(log(log(log(1/x)))) - 1)</pre>
<pre class="output">
f (generic function with 1 method)</pre>

<p>A numeric attempt might be done along these lines:</p><pre class="sourceCode julia">hs = [10.0^(-i) for i in 6:16]
ys = [f(h) for h in hs]
[hs ys]</pre>
<pre class="output">
11×2 Array{Float64,2}:
 1.0e-6   6.14632e-7 
 1.0e-7   1.42981e-7 
 1.0e-8   3.43858e-8 
 1.0e-9   8.52992e-9 
 1.0e-10  2.17687e-9 
 1.0e-11  5.70097e-10
 1.0e-12  1.52866e-10
 1.0e-13  4.18839e-11
 1.0e-14  1.17057e-11
 1.0e-15  3.33197e-12
 1.0e-16  9.64641e-13</pre>

<p>With a values appearing to approach $0$. However, in fact these values will ultimately head  off to $\infty$:</p><pre class="sourceCode julia">limit(f(x), x, 0, dir="+")</pre>
<div class="well well-sm">
$$\infty$$</div>

<h3>Derivatives</h3><p>One <em>could</em> use limits to implement the definition of a derivative:</p><pre class="sourceCode julia">x, h = symbols("x,h")
f(x) = exp(x)*sin(x)
limit((f(x+h) - f(x)) / h, h, 0)</pre>
<div class="well well-sm">
$$e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )}$$</div>

<p>However, it would be pretty inefficient, as <code>SymPy</code> already does a great job with derivatives. The <code>diff</code> function implements this. The basic syntax is <code>diff&#40;ex, x&#41;</code> to find the first derivative in <code>x</code> of the expression in <code>ex</code>, or its generalization to $k$th derivatives with <code>diff&#40;ex, x, k&#41;</code>.</p><p>The same derivative computed above by a limit could be found with:</p><pre class="sourceCode julia">diff(f(x), x)</pre>
<div class="well well-sm">
$$e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )}$$</div>

<p>Similarly, we can compute other derivatives:</p><pre class="sourceCode julia">diff(x^x, x)</pre>
<div class="well well-sm">
$$x^{x} \left(\log{\left (x \right )} + 1\right)$$</div>

<p>Or</p><pre class="sourceCode julia">diff(exp(-x^2), x, 2)</pre>
<div class="well well-sm">
$$2 \left(2 x^{2} - 1\right) e^{- x^{2}}$$</div>

<p>As an alternate to specifying the number of derivatives, multiple variables can be passed to <code>diff</code>:</p><pre class="sourceCode julia">diff(exp(-x^2), x, x, x)     # same as diff(..., x, 3)</pre>
<div class="well well-sm">
$$4 x \left(- 2 x^{2} + 3\right) e^{- x^{2}}$$</div>

<p>This could include variables besides <code>x</code>.</p><p>The output is a simple expression, so <code>diff</code> can be composed with other functions, such as <code>solve</code>. For example, here we find the critical points where the derivative is $0$ of some rational function:</p><pre class="sourceCode julia">f(x) = (12x^2 - 1) / (x^3)
diff(f(x), x) |> solve</pre>
<div class="well well-sm">
\begin{bmatrix}- \frac{1}{2}\\\frac{1}{2}\end{bmatrix}</div>

<h4>Operator version</h4><p><code>SymPy</code> provides an "operator" version of <code>diff</code> for univariate functions for convenience (<code>diff&#40;f::Function,k&#61;1&#41;&#61;diff&#40;f&#40;x&#41;,x,k&#41;</code>):</p><pre class="sourceCode julia">f(x) = exp(x)*cos(x)
diff(f, 2)</pre>
<div class="well well-sm">
$$- 2 e^{x} \sin{\left (x \right )}$$</div>

<h4>Partial derivatives</h4><p>The <code>diff</code> function makes finding partial derivatives as easy as specifying the variable to differentiate in. This  example computes the mixed partials of an expression in <code>x</code> and <code>y</code>:</p><pre class="sourceCode julia">x,y = symbols("x,y")
ex = x^2*cos(y)
Sym[diff(ex,v1, v2) for v1 in [x,y], v2 in [x,y]]</pre>
<div class="well well-sm">
\begin{bmatrix}2 \cos{\left (y \right )}&- 2 x \sin{\left (y \right )}\\- 2 x \sin{\left (y \right )}&- x^{2} \cos{\left (y \right )}\end{bmatrix}</div>

<p>The extra <code>Sym</code>, of the form <code>T&#91;&#93;</code>, helps <code>Julia</code> resolve the type of the output.</p><h4>Unevaluated derivatives</h4><p>The <code>Derivative</code> function provides unevaluated derivatives, useful with differential equations and the output for unknown functions. Here is an example:</p><pre class="sourceCode julia">ex = Derivative(exp(x*y), x, y, 2)</pre>
<div class="well well-sm">
$$\frac{\partial^{3}}{\partial x\partial y^{2}}  e^{x y}$$</div>

<p>(The <code>y,2</code> is a replacement for <code>y,y</code> which makes higher order terms easier to type.) These expressions are evaluated with <code>doit</code>:</p><pre class="sourceCode julia">doit(ex)</pre>
<div class="well well-sm">
$$x \left(x y + 2\right) e^{x y}$$</div>

<h4>Implicit derivatives</h4><p>SymPy can be used to find derivatives of implicitly defined functions. For example, the task of finding $dy/dx$ for the equation:</p>$$~
y^4 - x^4 -y^2 + 2x^2 = 0
~$$
<p>As with the mathematical solution, the key is to treat one of the variables as depending on the other. In this case, we think of $y$ locally as a function of $x$. SymPy allows us to create symbolic functions, and we will use one to substitute in for <code>y</code>.</p><p>In SymPy, symbolic functions use the class name  "Function", but in <code>SymPy</code> we use <code>SymFunction</code> to avoid a name collision with one of <code>Julia</code>'s primary types. The constructor can be used as <code>SymFunction&#40;:F&#41;</code>:</p><pre class="sourceCode julia">F, G = SymFunction("F"), SymFunction("G")</pre>
<pre class="output">
(,)</pre>

<p>We can call these functions, but we get a function expression:</p><pre class="sourceCode julia">F(x)</pre>
<div class="well well-sm">
$$F{\left (x \right )}$$</div>

<p>SymPy can differentiate symbolically, again with <code>diff</code>:</p><pre class="sourceCode julia">diff(F(x))</pre>
<div class="well well-sm">
$$\frac{d}{d x} F{\left (x \right )}$$</div>

<p>Of for symbolic functions the more natural <code>F&#39;&#40;x&#41;</code>.</p><p>To get back to our problem, we have our expression:</p><pre class="sourceCode julia">x,y = symbols("x, y")
ex = y^4 - x^4 - y^2 + 2x^2</pre>
<div class="well well-sm">
$$- x^{4} + 2 x^{2} + y^{4} - y^{2}$$</div>

<p>Now we substitute:</p><pre class="sourceCode julia">ex1 = ex(y=>F(x))</pre>
<div class="well well-sm">
$$- x^{4} + 2 x^{2} + F^{4}{\left (x \right )} - F^{2}{\left (x \right )}$$</div>

<p>We want to differentiate "both" sides. As the right side is just $0$, there isn't anything to do here, but mentally keep track. As for the left we have:</p><pre class="sourceCode julia">ex2 = diff(ex1, x)</pre>
<div class="well well-sm">
$$- 4 x^{3} + 4 x + 4 F^{3}{\left (x \right )} \frac{d}{d x} F{\left (x \right )} - 2 F{\left (x \right )} \frac{d}{d x} F{\left (x \right )}$$</div>

<p>Now we collect terms and solve in terms of $F'(x)$</p><pre class="sourceCode julia">ex3 = solve(ex2, F'(x))[1]</pre>
<div class="well well-sm">
$$\frac{2 x \left(x^{2} - 1\right)}{\left(2 F^{2}{\left (x \right )} - 1\right) F{\left (x \right )}}$$</div>

<p>Finally, we substitute back into the solution for $F(x)$:</p><pre class="sourceCode julia">ex4 = ex3(F(x) => y)</pre>
<div class="well well-sm">
$$\frac{2 x \left(x^{2} - 1\right)}{y \left(2 y^{2} - 1\right)}$$</div>

<h6>Example: A Norman Window</h6><p>A classic calculus problem is to maximize the area of a <a href="http://en.wiktionary.org/wiki/Norman_window">Norman window</a> (in the shape of a rectangle with a half circle atop) when the perimeter is fixed to be $P \geq 0$.</p><p>Label the rectangle with $w$ and $h$ for width and height and then the half circle has radius $r=w/2$. With this, we can see that the area is $wh+(1/2)\pi r^2$ and the perimeter is $w + 2h + \pi r$. This gives:</p><pre class="sourceCode julia">w, h, P = symbols("w, h, P", nonnegative=true)
r = w/2
A = w*h + 1//2 * (pi * r^2)
p = w + 2h + pi*r</pre>
<div class="well well-sm">
$$2 h + w + \frac{\pi w}{2}$$</div>

<p>(There is a subtlety above, as m <code>1//2*pi*r^2</code> will lose exactness, as the products will be done left to right, and <code>1//2*pi</code> will be converted to an approximate floating point value before multiplying <code>r^2</code>, as such we rewrite the terms. It may be easier to use <code>PI</code> instead of <code>pi</code>.)</p><p>We want to solve for <code>h</code> from when <code>p&#61;P</code> (our fixed value) and substitute back into <code>A</code>. We solve <code>P-p&#61;&#61;0</code>:</p><pre class="sourceCode julia">h0 =  solve(P-p, h)[1]
A1 = A(h => h0)</pre>
<div class="well well-sm">
$$\frac{\pi w^{2}}{8} + w \left(\frac{P}{2} - \frac{\pi w}{4} - \frac{w}{2}\right)$$</div>

<p>Now we note this is a parabola in <code>w</code>, so any maximum will be an endpoint or the vertex, provided the leading term is negative.  The leading term can be found through:</p><pre class="sourceCode julia">coeffs(Poly(A1, w))</pre>
<pre class="output">
2-element Array{Any,1}:
 -1/2 - pi/8
         P/2</pre>

<p>Or without using the <code>Poly</code> methods, we could do this:</p><pre class="sourceCode julia">coeff(collect(expand(A1), w), w^2)</pre>
<div class="well well-sm">
$$- \frac{1}{2} - \frac{\pi}{8}$$</div>

<p>Either way, the leading coefficient, $-1/2 - \pi/8$, is negative, so the maximum can only happen at an endpoint or the vertex of the parabola. Now we check that when $w=0$ (the left endpoint) the area is $0$:</p><pre class="sourceCode julia">A1(w => 0)</pre>
<div class="well well-sm">
$$0$$</div>

<p>The other endpoint is when $h=0$, or</p><pre class="sourceCode julia">b = solve(subs(P-p, h, 0), w)[1]</pre>
<div class="well well-sm">
$$\frac{2 P}{2 + \pi}$$</div>

<p>We will need to check the area at <code>b</code> and at the vertex.</p><p>To find the vertex, we can use calculus – it will be when the derivative in <code>w</code> is $0$:</p><pre class="sourceCode julia">c = solve(diff(A1, w), w)[1]</pre>
<div class="well well-sm">
$$\frac{2 P}{\pi + 4}$$</div>

<p>The answer will be the larger of <code>A1</code> at <code>b</code> or <code>c</code>:</p><pre class="sourceCode julia">atb = A1(w => b)
atc = A1(w => c)</pre>
<div class="well well-sm">
$$\frac{\pi P^{2}}{2 \left(\pi + 4\right)^{2}} + \frac{2 P}{\pi + 4} \left(- \frac{\pi P}{2 \pi + 8} - \frac{P}{\pi + 4} + \frac{P}{2}\right)$$</div>

<p>A simple comparison isn't revealing:</p><pre class="sourceCode julia">atc - atb</pre>
<div class="well well-sm">
$$- \frac{\pi P^{2}}{2 \left(2 + \pi\right)^{2}} + \frac{\pi P^{2}}{2 \left(\pi + 4\right)^{2}} - \frac{2 P}{2 + \pi} \left(- \frac{\pi P}{4 + 2 \pi} - \frac{P}{2 + \pi} + \frac{P}{2}\right) + \frac{2 P}{\pi + 4} \left(- \frac{\pi P}{2 \pi + 8} - \frac{P}{\pi + 4} + \frac{P}{2}\right)$$</div>

<p>But after simplifying, we can see that this expression is positive if $P$ is:</p><pre class="sourceCode julia">simplify(atc - atb)</pre>
<div class="well well-sm">
$$\frac{2 P^{2}}{16 + \pi^{3} + 20 \pi + 8 \pi^{2}}$$</div>

<p>With this observation, we conclude the maximum area happens at <code>c</code> with area <code>atc</code>.</p><h3>Integrals</h3><p>Integration is implemented in SymPy through the <code>integrate</code> function. There are two basic calls:  <code>integrate&#40;f&#40;x&#41;, x&#41;</code> will find the indefinite integral ($\int f(x) dx$) and when endpoints are specified through <code>integrate&#40;f&#40;x&#41;, &#40;x, a, b&#41;&#41;</code> the definite integral will be found ($\int_a^b f(x) dx$). The special form <code>integrate&#40;ex, x, a, b&#41;</code> can be used for single integrals, but the specification through a tuple is needed for multiple integrals.</p><p>Basic integrals are implemented:</p><pre class="sourceCode julia">integrate(x^3, x)</pre>
<div class="well well-sm">
$$\frac{x^{4}}{4}$$</div>

<p>Or in more generality:</p><pre class="sourceCode julia">n = symbols("n", real=true)
ex = integrate(x^n, x)</pre>
<div class="well well-sm">
$$\begin{cases} \log{\left (x \right )} & \text{for}\: n = -1 \\\frac{x^{n + 1}}{n + 1} & \text{otherwise} \end{cases}$$</div>

<p>The output here is a <em>piecewise function</em>, performing a substitution will choose a branch in this case:</p><pre class="sourceCode julia">ex(n => 3)</pre>
<div class="well well-sm">
$$\frac{x^{4}}{4}$$</div>

<p>Definite integrals are just as easy. Here is Archimedes' answer:</p><pre class="sourceCode julia">integrate(x^2, (x, 0, 1))</pre>
<div class="well well-sm">
$$\frac{1}{3}$$</div>

<p>Tedious problems, such as those needing multiple integration-by-parts steps can be done easily:</p><pre class="sourceCode julia">integrate(x^5*sin(x), x)</pre>
<div class="well well-sm">
$$- x^{5} \cos{\left (x \right )} + 5 x^{4} \sin{\left (x \right )} + 20 x^{3} \cos{\left (x \right )} - 60 x^{2} \sin{\left (x \right )} - 120 x \cos{\left (x \right )} + 120 \sin{\left (x \right )}$$</div>

<p>The SymPy tutorial says:</p><blockquote>
<p>"<code>integrate</code> uses powerful algorithms that are always improving to compute both definite and indefinite integrals, including heuristic pattern matching type algorithms, a partial implementation of the Risch algorithm, and an algorithm using Meijer G-functions that is useful for computing integrals in terms of special functions, especially definite integrals."</p>
</blockquote><p>The tutorial gives the following example:</p><pre class="sourceCode julia">f(x) = (x^4 + x^2 * exp(x) - x^2 - 2x*exp(x) - 2x - exp(x)) * exp(x) / ( (x-1)^2 * (x+1)^2 * (exp(x) + 1) )
f(x)</pre>
<div class="well well-sm">
$$\frac{\left(x^{4} + x^{2} e^{x} - x^{2} - 2 x e^{x} - 2 x - e^{x}\right) e^{x}}{\left(x - 1\right)^{2} \left(x + 1\right)^{2} \left(e^{x} + 1\right)}$$</div>

<p>With indefinite integral:</p><pre class="sourceCode julia">integrate(f(x), x)</pre>
<div class="well well-sm">
$$\log{\left (e^{x} + 1 \right )} + \frac{e^{x}}{x^{2} - 1}$$</div>

<h4>Multiple integrals</h4><p>The <code>integrate</code> function uses a tuple, <code>&#40;var, a, b&#41;</code>, to specify the limits of a definite integral. This syntax lends itself readily to multiple integration.</p><p>For example, the following computes the integral of $xy$ over the unit square:</p><pre class="sourceCode julia">x, y = symbols("x,y")
integrate(x*y, (y, 0, 1), (x, 0, 1))</pre>
<div class="well well-sm">
$$\frac{1}{4}$$</div>

<p>The innermost terms can depend on outer ones. For example, the following integrates $x^2y$ over the upper half of the unit circle:</p><pre class="sourceCode julia">integrate(x^2*y, (y, 0, sqrt(1 - x^2)), (x, -1, 1))</pre>
<div class="well well-sm">
$$\frac{2}{15}$$</div>

<h4>Unevaluated integrals</h4><p>The <code>Integral</code> function can stage unevaluated integrals that will be evaluated by calling <code>doit</code>. It is also used when the output is unknown. This example comes from the tutorial:</p><pre class="sourceCode julia">integ = Integral(sin(x^2), x)</pre>
<div class="well well-sm">
$$\int \sin{\left (x^{2} \right )}\, dx$$</div>

<pre class="sourceCode julia">doit(integ)</pre>
<div class="well well-sm">
$$\frac{3 \sqrt{2} \sqrt{\pi} S\left(\frac{\sqrt{2} x}{\sqrt{\pi}}\right)}{8 \Gamma{\left(\frac{7}{4} \right)}} \Gamma{\left(\frac{3}{4} \right)}$$</div>

<h4>Operator version</h4><p>For convenience, for univariate functions there is a convenience wrapper so that the operator styles – <code>integrate&#40;f&#41;</code> and <code>integrate&#40;f, a, b&#41;</code> – will perform the integrations.</p><pre class="sourceCode julia">f(x) = exp(x) * cos(x)
integrate(f)</pre>
<div class="well well-sm">
$$\frac{e^{x}}{2} \sin{\left (x \right )} + \frac{e^{x}}{2} \cos{\left (x \right )}$$</div>

<p>Or</p><pre class="sourceCode julia">integrate(sin, 0, pi)</pre>
<div class="well well-sm">
$$2$$</div>

<h3>Taylor series</h3><p>The <code>series</code> function can compute series expansions around a point to a specified order. For example, the following command finds 4 terms of the series expansion of <code>exp&#40;sin&#40;x&#41;&#41;</code> in <code>x</code> about $c=0$:</p><pre class="sourceCode julia">s1 = series(exp(sin(x)), x, 0, 4)</pre>
<div class="well well-sm">
$$1 + x + \frac{x^{2}}{2} + \mathcal{O}\left(x^{4}\right)$$</div>

<p>The coefficients are from the Taylor expansion ($a_i=f^{i}(c)/i!$). The <a href="http://en.wikipedia.org/wiki/Big_O_notation">big "O"</a> term indicates that any other power is no bigger than a constant times $x^4$.</p><p>Consider what happens when we multiply series of different orders:</p><pre class="sourceCode julia">s2 = series(cos(exp(x)), x, 0, 6)</pre>
<div class="well well-sm">
$$\cos{\left (1 \right )} - x \sin{\left (1 \right )} + x^{2} \left(- \frac{1}{2} \sin{\left (1 \right )} - \frac{1}{2} \cos{\left (1 \right )}\right) - \frac{x^{3}}{2} \cos{\left (1 \right )} + x^{4} \left(- \frac{1}{4} \cos{\left (1 \right )} + \frac{5}{24} \sin{\left (1 \right )}\right) + x^{5} \left(- \frac{1}{24} \cos{\left (1 \right )} + \frac{23}{120} \sin{\left (1 \right )}\right) + \mathcal{O}\left(x^{6}\right)$$</div>

<pre class="sourceCode julia">simplify(s1 * s2)</pre>
<div class="well well-sm">
$$\cos{\left (1 \right )} + \sqrt{2} x \cos{\left (\frac{\pi}{4} + 1 \right )} - \frac{3 x^{2}}{2} \sin{\left (1 \right )} - \sqrt{2} x^{3} \sin{\left (\frac{\pi}{4} + 1 \right )} + \mathcal{O}\left(x^{4}\right)$$</div>

<p>The big "O" term is $x^4$, as smaller order terms in <code>s2</code> are covered in this term. The big "O" notation is sometimes not desired, in which case the <code>removeO</code> function can be employed:</p><pre class="sourceCode julia">removeO(s1)</pre>
<div class="well well-sm">
$$\frac{x^{2}}{2} + x + 1$$</div>

<h3>Sums</h3><p><code>SymPy</code> can do sums, including some infinite ones. The <code>summation</code> function performs this task. For example, we have</p><pre class="sourceCode julia">i, n = symbols("i, n")
summation(i^2, (i, 1, n))</pre>
<div class="well well-sm">
$$\frac{n^{3}}{3} + \frac{n^{2}}{2} + \frac{n}{6}$$</div>

<p>Like <code>Integrate</code> and <code>Derivative</code>, there is also a <code>Sum</code> function to stage the task until the <code>doit</code> function is called to initiate the sum.</p><p>Some famous sums can be computed:</p><pre class="sourceCode julia">sn = Sum(1/i^2, (i, 1, n))
doit(sn)</pre>
<div class="well well-sm">
$$\operatorname{harmonic}{\left (n,2 \right )}$$</div>

<p>And from this a limit is available:</p><pre class="sourceCode julia">limit(doit(sn), n, oo)</pre>
<div class="well well-sm">
$$\frac{\pi^{2}}{6}$$</div>

<p>This would have also been possible through <code>summation&#40;1/i^2, &#40;i, 1, oo&#41;&#41;</code>.</p><h3>Vector-valued functions</h3><p>Julia makes constructing a vector of symbolic objects easy:</p><pre class="sourceCode julia">x,y = symbols("x,y")
v = [1,2,x]
w = [1,y,3]</pre>
<div class="well well-sm">
\begin{bmatrix}1\\y\\3\end{bmatrix}</div>

<p>The generic definitions of vector operations will work as expected with symbolic objects:</p><pre class="sourceCode julia">dot(v,w)</pre>
<div class="well well-sm">
$$2 y + 3 \overline{x} + 1$$</div>

<p>Or</p><pre class="sourceCode julia">cross(v,w)</pre>
<div class="well well-sm">
\begin{bmatrix}- x y + 6\\x - 3\\y - 2\end{bmatrix}</div>

<p>Finding gradients can be done using a comprehension.</p><pre class="sourceCode julia">ex = x^2*y - x*y^2
Sym[diff(ex,var) for var in [x,y]]</pre>
<div class="well well-sm">
\begin{bmatrix}2 x y - y^{2}\\x^{2} - 2 x y\end{bmatrix}</div>

<p>The mixed partials is similarly done by passing two variables to differentiate in to <code>diff</code>:</p><pre class="sourceCode julia">Sym[diff(ex, v1, v2) for v1 in [x,y], v2 in [x,y]]</pre>
<div class="well well-sm">
\begin{bmatrix}2 y&2 \left(x - y\right)\\2 \left(x - y\right)&- 2 x\end{bmatrix}</div>

<p>For this task, SymPy provides the <code>hessian</code> function:</p><pre class="sourceCode julia">hessian(ex)</pre>
<div class="well well-sm">
\begin{bmatrix}2 y&2 x - 2 y\\2 x - 2 y&- 2 x\end{bmatrix}</div>

<p>(When there are symbolic parameters, the free variables are specified as a vector, as in <code>hessian&#40;ex, vars&#41;</code>.)</p><h2>Matrices</h2><p>SymPy has a special class to work with matrices, as does <code>Julia</code>. With <code>SymPy</code>, matrices are just <code>Julia</code>n matrices with symbolic entries. The conversion to matrices that SymPy knows about is  handled in the background.</p><p>Constructing matrices then follows <code>Julia</code>'s conventions:</p><pre class="sourceCode julia">x,y = symbols("x,y")
M = [1 x; x 1]</pre>
<div class="well well-sm">
\begin{bmatrix}1&x\\x&1\end{bmatrix}</div>

<p>As much as possible, generic <code>Julia</code> functions are utilized:</p><pre class="sourceCode julia">diagm(ones(Sym, 5))
M^2
det(M)</pre>
<div class="well well-sm">
$$- x^{2} + 1$$</div>

<p>Occasionally, the SymPy method has more content:</p><pre class="sourceCode julia">eigvecs(M)</pre>
<pre class="output">
2-element Array{Any,1}:
 
⎡-1⎤
⎢  ⎥
⎣1 ⎦
 
⎡1⎤
⎢ ⎥
⎣1⎦   </pre>

<p>As compared to SymPy's <code>:egienvects</code> which yields:</p><pre class="sourceCode julia">M[:eigenvects]()</pre>
<pre class="output">
2-element Array{Any,1}:
 (-x + 1,1,Any[
⎡-1⎤
⎢  ⎥
⎣1 ⎦])
 (x + 1,1,Any[
⎡1⎤
⎢ ⎥
⎣1⎦])    </pre>

<p>This example from the tutorial shows the <code>nullspace</code> function:</p><pre class="sourceCode julia">M = Sym[1 2 3 0 0; 4 10 0 0 1]
vs = nullspace(M)</pre>
<pre class="output">
3-element Array{Any,1}:
 
⎡-15⎤
⎢   ⎥
⎢ 6 ⎥
⎢   ⎥
⎢ 1 ⎥
⎢   ⎥
⎢ 0 ⎥
⎢   ⎥
⎣ 0 ⎦         
 
⎡0⎤
⎢ ⎥
⎢0⎥
⎢ ⎥
⎢0⎥
⎢ ⎥
⎢1⎥
⎢ ⎥
⎣0⎦                           
 
⎡ 1  ⎤
⎢    ⎥
⎢-1/2⎥
⎢    ⎥
⎢ 0  ⎥
⎢    ⎥
⎢ 0  ⎥
⎢    ⎥
⎣ 1  ⎦</pre>

<p>And this shows that they are indeed in the null space of <code>M</code>:</p><pre class="sourceCode julia">[M*vs[i] for i in 1:3]</pre>
<pre class="output">
3-element Array{Array{SymPy.Sym,1},1}:
 
⎡0⎤
⎢ ⎥
⎣0⎦
 
⎡0⎤
⎢ ⎥
⎣0⎦
 
⎡0⎤
⎢ ⎥
⎣0⎦</pre>

<p>Symbolic expressions can be included in the matrices:</p><pre class="sourceCode julia">M = [1 x; x 1]
P, D = diagonalize(M)  # M = PDP^-1
D, M - P*D*inv(P)</pre>
<pre class="output">
(

⎡-x + 1    0  ⎤
⎢             ⎥
⎣  0     x + 1⎦,


⎡0  0⎤
⎢    ⎥
⎣0  0⎦)</pre>

<h2>Differential equations</h2><p>SymPy has facilities for solving ordinary differential <a href="http://docs.sympy.org/latest/modules/solvers/ode.html">equations</a>. The key is to create a symbolic function expression using  <code>SymFunction</code>. Again, this may be done through:</p><pre class="sourceCode julia">F = SymFunction("F")</pre>
<pre class="output">
F</pre>

<p>With this, we can  construct a  differential equation. Following the SymPy tutorial, we solve $f''(x) - 2f'(x) + f(x) = \sin(x)$:</p><pre class="sourceCode julia">diffeq = Eq(diff(F(x), x, 2) - 2*diff(F(x)) + F(x), sin(x))</pre>
<div class="well well-sm">
$$F{\left (x \right )} - 2 \frac{d}{d x} F{\left (x \right )} + \frac{d^{2}}{d x^{2}}  F{\left (x \right )} = \sin{\left (x \right )}$$</div>

<p>With this, we just need the <code>dsolve</code> function. This is called as <code>dsolve&#40;eq&#41;</code>:</p><pre class="sourceCode julia">ex = dsolve(diffeq)</pre>
<div class="well well-sm">
$$F{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{x} + \frac{1}{2} \cos{\left (x \right )}$$</div>

<p>The <code>dsolve</code> function in SymPy has an extensive list of named arguments to control the underlying algorithm. These can be passed through with the appropriate keyword arguments.</p><p>More clearly, the <code>SymFunction</code> objects have the <code>&#39;</code> method defined to find a derivative, so the above could also have been:</p><pre class="sourceCode julia">diffeq = F''(x) - 2F'(x) + F(x) - sin(x)
dsolve(diffeq)</pre>
<div class="well well-sm">
$$F{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{x} + \frac{1}{2} \cos{\left (x \right )}$$</div>

<p>This solution has two constants, $C_1$ and $C_2$, that would be found from initial conditions. Say we know $F(0)=0$ and $F'(0)=1$, can we find the constants? To work with the returned expression, it is most convenient to get just the right hand side. The <code>rhs</code> function will return the right-hand side of a relation:</p><pre class="sourceCode julia">ex1 = rhs(ex)</pre>
<div class="well well-sm">
$$\left(C_{1} + C_{2} x\right) e^{x} + \frac{1}{2} \cos{\left (x \right )}$$</div>

<p>(The <a href="http://docs.sympy.org/dev/modules/core.html#sympy.core.basic.Basic.args">args</a> function also can be used to break up the expression into parts.)</p><p>With this, we can solve for <code>C1</code> through substituting in $0$ for $x$:</p><pre class="sourceCode julia">solve(ex1(x => 0), Sym("C1"))</pre>
<div class="well well-sm">
\begin{bmatrix}- \frac{1}{2}\end{bmatrix}</div>

<p>We see that $C1=-1/2$, which we substitute in:</p><pre class="sourceCode julia">ex2 = ex1(Sym("C1") => -1//2)</pre>
<div class="well well-sm">
$$\left(C_{2} x - \frac{1}{2}\right) e^{x} + \frac{1}{2} \cos{\left (x \right )}$$</div>

<p>We know that $F'(0)=1$ now, so we solve for <code>C2</code> through</p><pre class="sourceCode julia">solve( subs(diff(ex2, x), x, 0) - 1, Sym("C2") )</pre>
<div class="well well-sm">
\begin{bmatrix}\frac{3}{2}\end{bmatrix}</div>

<p>This gives <code>C2&#61;3/2</code>. Again we substitute in to get our answer:</p><pre class="sourceCode julia">ex3 = ex2(Sym("C2") => 3//2)</pre>
<div class="well well-sm">
$$\left(\frac{3 x}{2} - \frac{1}{2}\right) e^{x} + \frac{1}{2} \cos{\left (x \right )}$$</div>

<h6>Example</h6><p>We do one more example, this one borrowed from <a href="http://nbviewer.ipython.org/github/garth-wells/IA-maths-Ipython/blob/master/notebooks/Lecture1.ipynb">here</a>.</p><blockquote>
<p>Find the variation of speed with time of a parachutist subject to a drag force of $k\cdot v^2$.</p>
</blockquote><p>The equation is</p>$$~
\frac{m}{k} \frac{dv}{dt} = \alpha^2 - v^2.
~$$
<p>We proceed through:</p><pre class="sourceCode julia">t, m,k,alpha = symbols("t,m,k,alpha")
v = SymFunction("v")
ex = Eq( (m/k)*v'(t), alpha^2 - v(t)^2 )</pre>
<div class="well well-sm">
$$\frac{m}{k} \frac{d}{d t} v{\left (t \right )} = \alpha^{2} - v^{2}{\left (t \right )}$$</div>

<p>We can "classify" this ODE with the method <code>classify_ode</code>. As this is not exported, we call it using indexing:</p><pre class="sourceCode julia">ex[:classify_ode]()</pre>
<pre class="output">
("separable","1st_power_series","lie_group","separable_Integral")</pre>

<p>It is linear, but not solvable. Proceeding with <code>dsolve</code> gives:</p><pre class="sourceCode julia">dsolve(ex)</pre>
<div class="well well-sm">
$$v{\left (t \right )} = - \frac{\alpha}{\tanh{\left (\frac{\alpha k}{m} \left(C_{1} - t\right) \right )}}$$</div>

<h3>Initial Value Problems</h3><p>Solving an initial value problem can be a bit tedious with <code>SymPy</code>. The first example shows the steps. This is because the <code>ics</code> argument for <code>dsolve</code> only works for a few types of equations. These do not include, by default, the familiar "book" examples, such as $y'(x) = a\cdot y(x)$.</p><p>To work around this, <code>SymPy.jl</code> extends the function <code>dsolve</code> to allow a specification of the initial conditions when solving.  The new ingredients are the independent variable (<code>x</code> in the examples) and tuples to specify each condition. The are conditions on the values of <code>u</code>, <code>u&#39;</code>', ....  To illustrate, we follow an example from <a href="https://reference.wolfram.com/language/tutorial/DSolveLinearBVPs.html">Wolfram</a>.</p><pre class="sourceCode julia">y = SymFunction("y")
a, x = symbols("a,x")
eqn = y'(x) - 3*x*y(x) - 1</pre>
<div class="well well-sm">
$$- 3 x y{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} - 1$$</div>

<p>We solve the initial value problem with $y(0) = 4$ as follows:</p><pre class="sourceCode julia">x0, y0 = 0, 4
out = dsolve(eqn, x, (y, x0, y0))</pre>
<div class="well well-sm">
$$y{\left (x \right )} = \left(\frac{\sqrt{6} \sqrt{\pi}}{6} \operatorname{erf}{\left (\frac{\sqrt{6} x}{2} \right )} + 4\right) e^{\frac{3 x^{2}}{2}}$$</div>

<p>Verifying this requires combining some operations:</p><pre class="sourceCode julia">u = rhs(out)
diff(u, x) - 3*x*u - 1</pre>
<div class="well well-sm">
$$0$$</div>

<p>To solve with a general initial condition is similar:</p><pre class="sourceCode julia">x0, y0 = 0, a
out = dsolve(eqn, x, (y, x0, y0))</pre>
<div class="well well-sm">
$$y{\left (x \right )} = \left(a + \frac{\sqrt{6} \sqrt{\pi}}{6} \operatorname{erf}{\left (\frac{\sqrt{6} x}{2} \right )}\right) e^{\frac{3 x^{2}}{2}}$$</div>

<p>To plot this over a range of values for <code>a</code> we have:</p><pre class="sourceCode julia">as = -2:0.6:2
ex = rhs(out)
p = plot(ex(a=>as[1]), -1.8, 1.8, ylims=(-4, 4))
for i in as[2:end]
  plot!(p, ex(a=>i), -1.8, 1.8, ylims=(-4, 4))
end
p  </pre>
    <div id="f1191879-e5c4-40ec-afce-89415b9268e0" style="width:600px;height:400px;"></div>
    <script>
    PLOT = document.getElementById('f1191879-e5c4-40ec-afce-89415b9268e0');
    Plotly.plot(PLOT, [{"yaxis":"y","y":[-351.24030439161226,-289.15883273572183,-238.99062515237088,-198.30588491902148,-165.19558215456863,-138.15505169961233,-115.99457813058983,-97.77050498700936,-82.73199879859288,-70.27978856888149,-59.93409265321788,-51.30961435109477,-44.09599166315742,-38.04246739003587,-32.945834072838046,-28.640927218027496,-24.99310695957327,-21.892295592582194,-19.248235851164118,-16.98670959481288,-15.046514127698973,-13.37703779033655,-11.936310824447174,-10.689434164069795,-9.607309531780036,-8.665610377207603,-7.843945824347734,-7.125179690008943,-6.494874409092808,-5.9408358236224075,-5.452739624515493,-5.021824058846146,-4.64063654850871,-4.302824278051386,-4.002960731862557,-3.736401696865121,-3.499165477028309,-3.2878330536049565,-3.099464719799844,-2.9315303597951434,-2.7818510606973548,-2.648550166578886,-2.530012225864803,-2.424848562386612,-2.331868428888213,-2.2500548895564405,-2.17854473324661,-2.1166118479332066,-2.0636535947775205,-2.0191798113552446,-1.982804151570118,-1.9542375375590049,-1.9332835589934894,-1.919835709801618,-1.9138764034421045,-1.9154777572932253,-1.9248041862408902,-1.9421168969540514,-1.9677804295068098,-2.0022714540134,-2.0461900991297766,-2.1002741693763216,-2.1654167025020015,-2.242687430473059,-2.3333588429403944,-2.4389377161593733,-2.561203170688048,-2.7022525670024296,-2.8645568510221513,-3.0510273360251325,-3.26509637198026,-3.510814931304073,-3.7929708611356086,-4.117232454238702,-4.490323121877141,-4.9202343742915,-5.416486107079062,-5.990445456905083,-6.655718359167908,-7.428631583627101,-8.328827662023933,-9.380001040711928,-10.610811364111882,-12.05601950808924,-13.757904472978812,-15.76803534941758,-18.149493385164295,-20.97966615642532,-24.353770897925347,-28.38930971004366,-33.231719011038585,-39.061553730152816,-46.10364934369434,-54.6388399768736,-65.0189892195224,-77.68632656139708,-93.1983960545533,-112.26034150377599,-135.76681020485302,-164.85650403790072],"showlegend":true,"name":"y1","type":"scatter","xaxis":"x","line":{"width":1,"dash":"solid","color":"rgba(0, 154, 250, 1.000)","shape":"linear"},"x":[-1.8,-1.7636363636363637,-1.7272727272727273,-1.690909090909091,-1.6545454545454545,-1.6181818181818182,-1.5818181818181818,-1.5454545454545454,-1.509090909090909,-1.4727272727272727,-1.4363636363636363,-1.4,-1.3636363636363635,-1.3272727272727274,-1.290909090909091,-1.2545454545454546,-1.2181818181818183,-1.1818181818181819,-1.1454545454545455,-1.1090909090909091,-1.0727272727272728,-1.0363636363636364,-1.0,-0.9636363636363636,-0.9272727272727272,-0.8909090909090909,-0.8545454545454545,-0.8181818181818182,-0.7818181818181819,-0.7454545454545455,-0.7090909090909091,-0.6727272727272727,-0.6363636363636364,-0.6,-0.5636363636363636,-0.5272727272727272,-0.4909090909090909,-0.45454545454545453,-0.41818181818181815,-0.38181818181818183,-0.34545454545454546,-0.3090909090909091,-0.2727272727272727,-0.23636363636363636,-0.2,-0.16363636363636364,-0.12727272727272726,-0.09090909090909091,-0.05454545454545454,-0.01818181818181818,0.01818181818181818,0.05454545454545454,0.09090909090909091,0.12727272727272726,0.16363636363636364,0.2,0.23636363636363636,0.2727272727272727,0.3090909090909091,0.34545454545454546,0.38181818181818183,0.41818181818181815,0.45454545454545453,0.4909090909090909,0.5272727272727272,0.5636363636363636,0.6,0.6363636363636364,0.6727272727272727,0.7090909090909091,0.7454545454545455,0.7818181818181819,0.8181818181818182,0.8545454545454545,0.8909090909090909,0.9272727272727272,0.9636363636363636,1.0,1.0363636363636364,1.0727272727272728,1.1090909090909091,1.1454545454545455,1.1818181818181819,1.2181818181818183,1.2545454545454546,1.290909090909091,1.3272727272727274,1.3636363636363635,1.4,1.4363636363636363,1.4727272727272727,1.509090909090909,1.5454545454545454,1.5818181818181818,1.6181818181818182,1.6545454545454545,1.690909090909091,1.7272727272727273,1.7636363636363637,1.8],"mode":"lines"},{"yaxis":"y","y":[-273.82578312718533,-225.4199862946356,-186.30298015394885,-154.58024277298526,-128.76329584717377,-107.67894556174213,-90.3995654144703,-76.1893818374038,-64.46296591928103,-54.753062431893476,-46.68558229872865,-39.96010656374175,-34.33464299022001,-29.613673273755843,-25.638753659499702,-22.28110246437655,-19.43573798942389,-17.016829549078082,-14.954000317382711,-13.189379006287357,-11.675242271000062,-10.372124367910882,-9.247297382244335,-8.273546123397466,-7.428177945869305,-6.692220352344891,-6.049769082559768,-5.487457107337259,-4.994021008033275,-4.559945994282007,-4.177174580434398,-3.8388669223659018,-3.5392031811820392,-3.273220160740471,-3.036675964659267,-2.8259376158942935,-2.6378875409031037,-2.4698455901889127,-2.3195038864234188,-2.1848722909564056,-2.0642326834907414,-1.9561005771660318,-1.8591928574419747,-1.7724006500924865,-1.6947665009609973,-1.6254651956066588,-1.5637876667893755,-1.5091275368942023,-1.4609699249270416,-1.4188822169164401,-1.3825065571313135,-1.351553867708526,-1.325799247954485,-1.3050786433443837,-1.2892867094923226,-1.2783758293660095,-1.2723562739467649,-1.2712975285312231,-1.2753308400939554,-1.2846530768067868,-1.2995320302910385,-1.3203133359998969,-1.3474292390859577,-1.3814094943478537,-1.4228947619695667,-1.4726529489560833,-1.531599053377133,-1.6008191996757588,-1.6815997145419064,-1.7754622919440386,-1.88420654263986,-2.0099615302445404,-2.1552482784639255,-2.3230557124507367,-2.516933097014429,-2.7411027883807693,-3.000598066406733,-3.3014320147022445,-3.6508049367422393,-4.057359726928189,-4.531497073498411,-5.085765506930521,-5.73534532060777,-6.498650537939861,-7.398079719327864,-8.460954936079236,-9.72069926888427,-11.218317483487908,-13.004263110572328,-15.140799355554424,-17.704992874050568,-20.792520850840962,-24.522526194088783,-29.043827260754078,-34.542883081652185,-41.25404025400222,-49.47275390851708,-59.572696505353946,-72.02796376376679,-87.44198277347377],"showlegend":true,"name":"y2","type":"scatter","xaxis":"x","line":{"width":1,"dash":"solid","color":"rgba(227, 111, 71, 1.000)","shape":"linear"},"x":[-1.8,-1.7636363636363637,-1.7272727272727273,-1.690909090909091,-1.6545454545454545,-1.6181818181818182,-1.5818181818181818,-1.5454545454545454,-1.509090909090909,-1.4727272727272727,-1.4363636363636363,-1.4,-1.3636363636363635,-1.3272727272727274,-1.290909090909091,-1.2545454545454546,-1.2181818181818183,-1.1818181818181819,-1.1454545454545455,-1.1090909090909091,-1.0727272727272728,-1.0363636363636364,-1.0,-0.9636363636363636,-0.9272727272727272,-0.8909090909090909,-0.8545454545454545,-0.8181818181818182,-0.7818181818181819,-0.7454545454545455,-0.7090909090909091,-0.6727272727272727,-0.6363636363636364,-0.6,-0.5636363636363636,-0.5272727272727272,-0.4909090909090909,-0.45454545454545453,-0.41818181818181815,-0.38181818181818183,-0.34545454545454546,-0.3090909090909091,-0.2727272727272727,-0.23636363636363636,-0.2,-0.16363636363636364,-0.12727272727272726,-0.09090909090909091,-0.05454545454545454,-0.01818181818181818,0.01818181818181818,0.05454545454545454,0.09090909090909091,0.12727272727272726,0.16363636363636364,0.2,0.23636363636363636,0.2727272727272727,0.3090909090909091,0.34545454545454546,0.38181818181818183,0.41818181818181815,0.45454545454545453,0.4909090909090909,0.5272727272727272,0.5636363636363636,0.6,0.6363636363636364,0.6727272727272727,0.7090909090909091,0.7454545454545455,0.7818181818181819,0.8181818181818182,0.8545454545454545,0.8909090909090909,0.9272727272727272,0.9636363636363636,1.0,1.0363636363636364,1.0727272727272728,1.1090909090909091,1.1454545454545455,1.1818181818181819,1.2181818181818183,1.2545454545454546,1.290909090909091,1.3272727272727274,1.3636363636363635,1.4,1.4363636363636363,1.4727272727272727,1.509090909090909,1.5454545454545454,1.5818181818181818,1.6181818181818182,1.6545454545454545,1.690909090909091,1.7272727272727273,1.7636363636363637,1.8],"mode":"lines"},{"yaxis":"y","y":[-196.41126186275835,-161.68113985354935,-133.61533515552682,-110.85460062694905,-92.33100953977892,-77.20283942387192,-64.80455269835079,-54.60825868779825,-46.193933039969174,-39.226336294905465,-33.43707194423942,-28.610598776388738,-24.5732943172826,-21.18487915747582,-18.33167324616136,-15.921277710725603,-13.878369019274519,-12.14136350557397,-10.659764783601304,-9.392048417761835,-8.303970414301151,-7.367210945485212,-6.558283940041497,-5.857658082725138,-5.249046359958575,-4.71883032748218,-4.255592340771803,-3.849734524665577,-3.4931676069737434,-3.1790561649416067,-2.901609536353305,-2.655909785885657,-2.437769813855368,-2.243616043429556,-2.0703911974559777,-1.9154735349234666,-1.7766096047778988,-1.651858126772869,-1.539543053046994,-1.4382142221176675,-1.3466143062841283,-1.2636509877531772,-1.1883734890191466,-1.1199527377983611,-1.0576645730337815,-1.000875501656877,-0.9490306003321416,-0.9016432258551977,-0.8582862550765629,-0.8185846224776357,-0.7822089626925092,-0.7488701978580472,-0.7183149369154805,-0.6903215768871496,-0.6646970155425409,-0.6412739014387938,-0.6199083616526395,-0.600478160108395,-0.5828812506811009,-0.5670346996001736,-0.5528739614523005,-0.5403525026234721,-0.5294417756699139,-0.5201315582226484,-0.5124306809987396,-0.5063681817527941,-0.5019949360662178,-0.49938583234908757,-0.49864257806166173,-0.49989724786294504,-0.5033167132994598,-0.5091081291850081,-0.5175256957922431,-0.5288789706627715,-0.5435430721517177,-0.5619712024700387,-0.5847100257344049,-0.6124185724994053,-0.6458915143165704,-0.6860878702292779,-0.7341664849728886,-0.7915299731491134,-0.8598792771036584,-0.9412815677904858,-1.0382549656769182,-1.1538745227408926,-1.291905152604245,-1.4569688105504965,-1.6547553232193097,-1.8922890010651967,-2.1782667370625624,-2.523487971529106,-2.9414030444832293,-3.448814544634566,-4.0667769437819805,-4.821753946607362,-5.7471117624808565,-6.885051506931918,-8.289117322680555,-10.02746150904681],"showlegend":true,"name":"y3","type":"scatter","xaxis":"x","line":{"width":1,"dash":"solid","color":"rgba(62, 164, 78, 1.000)","shape":"linear"},"x":[-1.8,-1.7636363636363637,-1.7272727272727273,-1.690909090909091,-1.6545454545454545,-1.6181818181818182,-1.5818181818181818,-1.5454545454545454,-1.509090909090909,-1.4727272727272727,-1.4363636363636363,-1.4,-1.3636363636363635,-1.3272727272727274,-1.290909090909091,-1.2545454545454546,-1.2181818181818183,-1.1818181818181819,-1.1454545454545455,-1.1090909090909091,-1.0727272727272728,-1.0363636363636364,-1.0,-0.9636363636363636,-0.9272727272727272,-0.8909090909090909,-0.8545454545454545,-0.8181818181818182,-0.7818181818181819,-0.7454545454545455,-0.7090909090909091,-0.6727272727272727,-0.6363636363636364,-0.6,-0.5636363636363636,-0.5272727272727272,-0.4909090909090909,-0.45454545454545453,-0.41818181818181815,-0.38181818181818183,-0.34545454545454546,-0.3090909090909091,-0.2727272727272727,-0.23636363636363636,-0.2,-0.16363636363636364,-0.12727272727272726,-0.09090909090909091,-0.05454545454545454,-0.01818181818181818,0.01818181818181818,0.05454545454545454,0.09090909090909091,0.12727272727272726,0.16363636363636364,0.2,0.23636363636363636,0.2727272727272727,0.3090909090909091,0.34545454545454546,0.38181818181818183,0.41818181818181815,0.45454545454545453,0.4909090909090909,0.5272727272727272,0.5636363636363636,0.6,0.6363636363636364,0.6727272727272727,0.7090909090909091,0.7454545454545455,0.7818181818181819,0.8181818181818182,0.8545454545454545,0.8909090909090909,0.9272727272727272,0.9636363636363636,1.0,1.0363636363636364,1.0727272727272728,1.1090909090909091,1.1454545454545455,1.1818181818181819,1.2181818181818183,1.2545454545454546,1.290909090909091,1.3272727272727274,1.3636363636363635,1.4,1.4363636363636363,1.4727272727272727,1.509090909090909,1.5454545454545454,1.5818181818181818,1.6181818181818182,1.6545454545454545,1.690909090909091,1.7272727272727273,1.7636363636363637,1.8],"mode":"lines"},{"yaxis":"y","y":[-118.99674059833141,-97.94229341246314,-80.92769015710479,-67.12895848091283,-55.89872323238406,-46.7267332860017,-39.209539982231284,-33.027135538192695,-27.92490016065732,-23.699610157917455,-20.188561589750186,-17.26109098903572,-14.811945644345187,-12.756085041195794,-11.024592832823014,-9.561452957074657,-8.321000049125141,-7.2658974620698595,-6.365529249819898,-5.5947178292363136,-4.93269855760224,-4.362297523059543,-3.8692704978386585,-3.4417700420528092,-3.0699147740478447,-2.7454403026194685,-2.4614155989838373,-2.2120119419938944,-1.9923142059142116,-1.7981663356012068,-1.6260444922722113,-1.4729526494054124,-1.336336446528697,-1.2140119261186406,-1.1041064302526884,-1.0050094539526393,-0.9153316686526936,-0.8338706633568254,-0.7595822196705693,-0.6915561532789294,-0.628995929077515,-0.5712013983403229,-0.5175541205963186,-0.4675048255042359,-0.4205626451065657,-0.3762858077070953,-0.3342735338749074,-0.2941589148161934,-0.25560258522608403,-0.21828702803883135,-0.18191136825370482,-0.14618652800756843,-0.11083062587647614,-0.07556451042991538,-0.04010732159275914,-0.004171973511578042,0.0325395506414858,0.07034120831443316,0.10956833873175344,0.1505836776064396,0.1937841073864375,0.2396083307529528,0.2885456877461297,0.3411463779025568,0.3980333999720878,0.45991658545049546,0.5276091812446972,0.6020475349775832,0.6843145584185829,0.7756677962181489,0.8775731160409401,0.9917452718745238,1.1201968868794394,1.265297771125194,1.429846952710994,1.6171603834406914,1.8311780149379238,2.0765948697034333,2.359021908109098,2.6851839864696334,3.0631641035526327,3.502705560632293,4.015586766400452,4.616087402358891,5.321569787974027,6.153205890597452,7.136888963675779,8.304379862386915,9.694752464133709,11.356221353424035,13.34845939992545,15.745544907782751,18.639720105122326,22.146198171484944,26.409329194088233,31.610532360787495,37.97853038355536,45.80259349149011,55.44972911840567,67.38705975538014],"showlegend":true,"name":"y4","type":"scatter","xaxis":"x","line":{"width":1,"dash":"solid","color":"rgba(195, 113, 210, 1.000)","shape":"linear"},"x":[-1.8,-1.7636363636363637,-1.7272727272727273,-1.690909090909091,-1.6545454545454545,-1.6181818181818182,-1.5818181818181818,-1.5454545454545454,-1.509090909090909,-1.4727272727272727,-1.4363636363636363,-1.4,-1.3636363636363635,-1.3272727272727274,-1.290909090909091,-1.2545454545454546,-1.2181818181818183,-1.1818181818181819,-1.1454545454545455,-1.1090909090909091,-1.0727272727272728,-1.0363636363636364,-1.0,-0.9636363636363636,-0.9272727272727272,-0.8909090909090909,-0.8545454545454545,-0.8181818181818182,-0.7818181818181819,-0.7454545454545455,-0.7090909090909091,-0.6727272727272727,-0.6363636363636364,-0.6,-0.5636363636363636,-0.5272727272727272,-0.4909090909090909,-0.45454545454545453,-0.41818181818181815,-0.38181818181818183,-0.34545454545454546,-0.3090909090909091,-0.2727272727272727,-0.23636363636363636,-0.2,-0.16363636363636364,-0.12727272727272726,-0.09090909090909091,-0.05454545454545454,-0.01818181818181818,0.01818181818181818,0.05454545454545454,0.09090909090909091,0.12727272727272726,0.16363636363636364,0.2,0.23636363636363636,0.2727272727272727,0.3090909090909091,0.34545454545454546,0.38181818181818183,0.41818181818181815,0.45454545454545453,0.4909090909090909,0.5272727272727272,0.5636363636363636,0.6,0.6363636363636364,0.6727272727272727,0.7090909090909091,0.7454545454545455,0.7818181818181819,0.8181818181818182,0.8545454545454545,0.8909090909090909,0.9272727272727272,0.9636363636363636,1.0,1.0363636363636364,1.0727272727272728,1.1090909090909091,1.1454545454545455,1.1818181818181819,1.2181818181818183,1.2545454545454546,1.290909090909091,1.3272727272727274,1.3636363636363635,1.4,1.4363636363636363,1.4727272727272727,1.509090909090909,1.5454545454545454,1.5818181818181818,1.6181818181818182,1.6545454545454545,1.690909090909091,1.7272727272727273,1.7636363636363637,1.8],"mode":"lines"},{"yaxis":"y","y":[-41.58221933390449,-34.203446971376934,-28.240045158682765,-23.403316334876624,-19.46643692498921,-16.250627148131496,-13.614527266111775,-11.446012388587143,-9.655867281345467,-8.172884020929446,-6.940051235260959,-5.911583201682703,-5.050596971407778,-4.327290924915772,-3.7175124194846707,-3.201628203423712,-2.763631078975766,-2.390431418565749,-2.0712937160384914,-1.7973872407107923,-1.5614267009033298,-1.3573841006338756,-1.1802570556358205,-1.025882001380481,-0.8907831881371149,-0.7720502777567569,-0.6672388571958723,-0.5742893593222119,-0.49146080485467974,-0.4172765062608069,-0.3504794481911177,-0.28999551292516806,-0.2349030792020263,-0.18440780880772564,-0.137821663049399,-0.09454537298181208,-0.05405373252748846,-0.015883199940781828,0.020378613705855416,0.05510191555980845,0.08862244812909807,0.1212481910725313,0.15326524782650952,0.18494308678988924,0.21653928282064994,0.24830388624268637,0.2804835325823267,0.3133253962228109,0.3470810846243947,0.38201056639997294,0.4183862261850994,0.45649714184291024,0.4966536851625281,0.5391925560273187,0.5844823723570225,0.6329299544156376,0.684987462935611,0.7411605767372612,0.8020179281446077,0.8682020548130527,0.9404421762251753,1.0195691641293774,1.1065331511621732,1.2024243140277617,1.3084974809429148,1.4262013526537847,1.5572132985556122,1.703480902304254,1.8672716948988273,2.0512328402992424,2.25846294538134,2.4925986729340557,2.7579194695511218,3.059474512913159,3.4032369775737052,3.7962919693514214,4.247066055610252,4.765608311906271,5.3639353305347655,6.056455843168544,6.860494692078154,7.796941094413699,8.891052809904563,10.173456372508266,11.681394541624972,13.460286303935794,15.565683079955802,18.065728535324325,21.044260251486726,24.604731707913263,28.87518553691346,34.0145777870946,40.220843254727875,47.74121088760445,56.88543533195844,68.04281866818235,81.70417252959157,98.49023848991213,119.18857555949188,144.80158101980706],"showlegend":true,"name":"y5","type":"scatter","xaxis":"x","line":{"width":1,"dash":"solid","color":"rgba(172, 142, 24, 1.000)","shape":"linear"},"x":[-1.8,-1.7636363636363637,-1.7272727272727273,-1.690909090909091,-1.6545454545454545,-1.6181818181818182,-1.5818181818181818,-1.5454545454545454,-1.509090909090909,-1.4727272727272727,-1.4363636363636363,-1.4,-1.3636363636363635,-1.3272727272727274,-1.290909090909091,-1.2545454545454546,-1.2181818181818183,-1.1818181818181819,-1.1454545454545455,-1.1090909090909091,-1.0727272727272728,-1.0363636363636364,-1.0,-0.9636363636363636,-0.9272727272727272,-0.8909090909090909,-0.8545454545454545,-0.8181818181818182,-0.7818181818181819,-0.7454545454545455,-0.7090909090909091,-0.6727272727272727,-0.6363636363636364,-0.6,-0.5636363636363636,-0.5272727272727272,-0.4909090909090909,-0.45454545454545453,-0.41818181818181815,-0.38181818181818183,-0.34545454545454546,-0.3090909090909091,-0.2727272727272727,-0.23636363636363636,-0.2,-0.16363636363636364,-0.12727272727272726,-0.09090909090909091,-0.05454545454545454,-0.01818181818181818,0.01818181818181818,0.05454545454545454,0.09090909090909091,0.12727272727272726,0.16363636363636364,0.2,0.23636363636363636,0.2727272727272727,0.3090909090909091,0.34545454545454546,0.38181818181818183,0.41818181818181815,0.45454545454545453,0.4909090909090909,0.5272727272727272,0.5636363636363636,0.6,0.6363636363636364,0.6727272727272727,0.7090909090909091,0.7454545454545455,0.7818181818181819,0.8181818181818182,0.8545454545454545,0.8909090909090909,0.9272727272727272,0.9636363636363636,1.0,1.0363636363636364,1.0727272727272728,1.1090909090909091,1.1454545454545455,1.1818181818181819,1.2181818181818183,1.2545454545454546,1.290909090909091,1.3272727272727274,1.3636363636363635,1.4,1.4363636363636363,1.4727272727272727,1.509090909090909,1.5454545454545454,1.5818181818181818,1.6181818181818182,1.6545454545454545,1.690909090909091,1.7272727272727273,1.7636363636363637,1.8],"mode":"lines"},{"yaxis":"y","y":[35.83230193052247,29.535399469709304,24.447599839739272,20.3223258111596,16.965849382405654,14.225478989738717,11.980485450007741,10.135110761018415,8.613165597966391,7.3538421160585665,6.3084591192282735,5.437924585670317,4.710751701529635,4.101503191364254,3.589567993853674,3.1581965502272347,2.7937378911736115,2.485034624938363,2.222941817742916,1.99994334781473,1.8098451557955821,1.6475293217917937,1.5087563865670186,1.3900060392918479,1.2883483977736159,1.2013397471059548,1.1269378845920934,1.0634332233494708,1.0093925962048524,0.9636133230795934,0.9250855958899764,0.8929616235550768,0.8665302881246447,0.8451963085031896,0.8284631041538906,0.8159187079890153,0.8072242035977168,0.802104263475262,0.8003394470822804,0.8017599843985466,0.8062408253357114,0.8136977804853858,0.8240846162493378,0.8373909990840147,0.8536412107478658,0.8728935801924682,0.895240599039561,0.9208097072618154,0.9497647544748735,0.9823081608387774,1.018683820623904,1.0591808116933892,1.1041379962015327,1.153949622484553,1.2090720663068044,1.2700318823428536,1.3374353752297363,1.4119799451600894,1.4944675175574622,1.585820432019666,1.6871002450639134,1.7995299975058023,1.924520614578217,2.063702250152967,2.2189615619137424,2.3924861198570744,2.5868174158665274,2.8049142696309253,3.050228831379072,3.3267978843803365,3.6393527747217402,3.9934520739935877,4.395642052222804,4.853651254701124,5.376627002436417,5.975423555262152,6.662954096282581,7.454621754109111,8.368848752960435,9.427727699867456,10.657825280603676,12.091176628195107,13.766518853408675,15.730825342657644,18.04121929527592,20.76736671727414,23.994477196235827,27.827077208261738,32.393768038839745,37.853242062402494,44.40191167390147,52.28361066640646,61.801966404333434,73.33622360372397,87.36154146982865,104.47510497557721,125.4298146756278,151.17788348833417,182.9274220005781,222.216102284234],"showlegend":true,"name":"y6","type":"scatter","xaxis":"x","line":{"width":1,"dash":"solid","color":"rgba(0, 170, 174, 1.000)","shape":"linear"},"x":[-1.8,-1.7636363636363637,-1.7272727272727273,-1.690909090909091,-1.6545454545454545,-1.6181818181818182,-1.5818181818181818,-1.5454545454545454,-1.509090909090909,-1.4727272727272727,-1.4363636363636363,-1.4,-1.3636363636363635,-1.3272727272727274,-1.290909090909091,-1.2545454545454546,-1.2181818181818183,-1.1818181818181819,-1.1454545454545455,-1.1090909090909091,-1.0727272727272728,-1.0363636363636364,-1.0,-0.9636363636363636,-0.9272727272727272,-0.8909090909090909,-0.8545454545454545,-0.8181818181818182,-0.7818181818181819,-0.7454545454545455,-0.7090909090909091,-0.6727272727272727,-0.6363636363636364,-0.6,-0.5636363636363636,-0.5272727272727272,-0.4909090909090909,-0.45454545454545453,-0.41818181818181815,-0.38181818181818183,-0.34545454545454546,-0.3090909090909091,-0.2727272727272727,-0.23636363636363636,-0.2,-0.16363636363636364,-0.12727272727272726,-0.09090909090909091,-0.05454545454545454,-0.01818181818181818,0.01818181818181818,0.05454545454545454,0.09090909090909091,0.12727272727272726,0.16363636363636364,0.2,0.23636363636363636,0.2727272727272727,0.3090909090909091,0.34545454545454546,0.38181818181818183,0.41818181818181815,0.45454545454545453,0.4909090909090909,0.5272727272727272,0.5636363636363636,0.6,0.6363636363636364,0.6727272727272727,0.7090909090909091,0.7454545454545455,0.7818181818181819,0.8181818181818182,0.8545454545454545,0.8909090909090909,0.9272727272727272,0.9636363636363636,1.0,1.0363636363636364,1.0727272727272728,1.1090909090909091,1.1454545454545455,1.1818181818181819,1.2181818181818183,1.2545454545454546,1.290909090909091,1.3272727272727274,1.3636363636363635,1.4,1.4363636363636363,1.4727272727272727,1.509090909090909,1.5454545454545454,1.5818181818181818,1.6181818181818182,1.6545454545454545,1.690909090909091,1.7272727272727273,1.7636363636363637,1.8],"mode":"lines"},{"yaxis":"y","y":[113.24682319494943,93.27424591079554,77.13524483816131,64.04796795719582,53.39813568980051,44.70158512760893,37.57549816612726,31.716233910623977,26.882198477278255,22.880568253046583,19.556969473717505,16.787432373023336,14.472100374467047,12.53029730764428,10.896648407192018,9.518021303878182,8.35110686132299,7.360500668442475,6.5171773515243245,5.797273936340252,5.181117012494494,4.652442744217463,4.197769828769857,3.8058940799641774,3.4674799836843464,3.1747297719686665,2.921114626380059,2.701155806021154,2.510245997264385,2.3445031524199935,2.2006506399710704,2.0759187600353215,1.9679636554513158,1.874800425814105,1.79474787135718,1.7263827889598427,1.668502139722922,1.6200917268913058,1.5803002804587052,1.5484180532372847,1.5238592025423248,1.5061473698982402,1.4949039846721661,1.48983891137814,1.4907431386750816,1.4974832741422501,1.5099976654967955,1.5282940183008198,1.5524484243253525,1.5826057552775818,1.6189814150627084,1.6618644815438681,1.7116223072405372,1.7687066889417873,1.8336617602565863,1.9071338102700695,1.9898832875238617,2.082799313582918,2.1869171069703164,2.303438809226279,2.4337583139026515,2.5794908308822273,2.7425080779942608,2.924980186278172,3.12942564288457,3.3587708870603636,3.616421533177443,3.9063476369575962,4.233185967859317,4.602362928461431,5.02024260406214,5.49430547505312,6.033364634894488,6.647827996489091,7.3500170272991285,8.154555141172882,9.07884213695491,10.14363519631195,11.373762175386105,12.798999556566367,14.455155869129198,16.385412161976515,18.64198489691279,21.288194312807022,24.401044048926867,28.074447130612484,32.42327131251585,37.58842588119915,43.743275826192765,51.10175241689173,59.92863781088949,70.55264354571833,83.383089553939,98.93123631984349,117.83764760769887,140.90739128297207,169.15545682166402,203.8655284867562,246.66626844166436,299.630623548661],"showlegend":true,"name":"y7","type":"scatter","xaxis":"x","line":{"width":1,"dash":"solid","color":"rgba(237, 94, 147, 1.000)","shape":"linear"},"x":[-1.8,-1.7636363636363637,-1.7272727272727273,-1.690909090909091,-1.6545454545454545,-1.6181818181818182,-1.5818181818181818,-1.5454545454545454,-1.509090909090909,-1.4727272727272727,-1.4363636363636363,-1.4,-1.3636363636363635,-1.3272727272727274,-1.290909090909091,-1.2545454545454546,-1.2181818181818183,-1.1818181818181819,-1.1454545454545455,-1.1090909090909091,-1.0727272727272728,-1.0363636363636364,-1.0,-0.9636363636363636,-0.9272727272727272,-0.8909090909090909,-0.8545454545454545,-0.8181818181818182,-0.7818181818181819,-0.7454545454545455,-0.7090909090909091,-0.6727272727272727,-0.6363636363636364,-0.6,-0.5636363636363636,-0.5272727272727272,-0.4909090909090909,-0.45454545454545453,-0.41818181818181815,-0.38181818181818183,-0.34545454545454546,-0.3090909090909091,-0.2727272727272727,-0.23636363636363636,-0.2,-0.16363636363636364,-0.12727272727272726,-0.09090909090909091,-0.05454545454545454,-0.01818181818181818,0.01818181818181818,0.05454545454545454,0.09090909090909091,0.12727272727272726,0.16363636363636364,0.2,0.23636363636363636,0.2727272727272727,0.3090909090909091,0.34545454545454546,0.38181818181818183,0.41818181818181815,0.45454545454545453,0.4909090909090909,0.5272727272727272,0.5636363636363636,0.6,0.6363636363636364,0.6727272727272727,0.7090909090909091,0.7454545454545455,0.7818181818181819,0.8181818181818182,0.8545454545454545,0.8909090909090909,0.9272727272727272,0.9636363636363636,1.0,1.0363636363636364,1.0727272727272728,1.1090909090909091,1.1454545454545455,1.1818181818181819,1.2181818181818183,1.2545454545454546,1.290909090909091,1.3272727272727274,1.3636363636363635,1.4,1.4363636363636363,1.4727272727272727,1.509090909090909,1.5454545454545454,1.5818181818181818,1.6181818181818182,1.6545454545454545,1.690909090909091,1.7272727272727273,1.7636363636363637,1.8],"mode":"lines"}], {"yaxis":{"type":"-","titlefont":{"size":15,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"title":"","tickfont":{"size":11,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"ticks":"inside","range":[-4.0,4.0],"tickmode":"array","showgrid":true,"tickvals":[-4.0,-2.0,0.0,2.0,4.0],"domain":[0.03762029746281716,0.9901574803149606],"ticktext":["-4","-2","0","2","4"],"tickangle":0,"zeroline":false,"linecolor":"rgba(0, 0, 0, 1.000)","tickcolor":"rgba(0, 0, 0, 1.000)","anchor":"x"},"annotations":[],"width":600,"plot_bgcolor":"rgba(255, 255, 255, 1.000)","showlegend":true,"legend":{"bgcolor":"rgba(255, 255, 255, 1.000)","y":1.0,"font":{"size":11,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"bordercolor":"rgba(0, 0, 0, 1.000)","x":1.0},"xaxis":{"type":"-","titlefont":{"size":15,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"title":"","tickfont":{"size":11,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"ticks":"inside","tickmode":"array","showgrid":true,"tickvals":[-1.0,0.0,1.0],"domain":[0.03619130941965587,0.9934383202099738],"ticktext":["-1","0","1"],"tickangle":0,"zeroline":false,"linecolor":"rgba(0, 0, 0, 1.000)","tickcolor":"rgba(0, 0, 0, 1.000)","anchor":"y"},"paper_bgcolor":"rgba(255, 255, 255, 1.000)","height":400,"margin":{"r":0,"l":0,"b":0,"t":20}});
    </script>
<p>The comment from the example is "This plots several integral curves of the equation for different values of $a$. The plot shows that the solutions have an inflection point if the parameter  lies between $-1$ and $1$ , while a global maximum or minimum arises for other values of $a$."</p><h5>Example</h5><p>We continue with another example from the Wolfram documentation, that of solving $y'' + 5y' + 6y=0$ with values prescribed for both $y$ and $y'$ at $x_0=0$.</p><pre class="sourceCode julia">y = SymFunction("y")
x = symbols("x")
eqn = y''(x) + 5y'(x) + 6y(x)</pre>
<div class="well well-sm">
$$6 y{\left (x \right )} + 5 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}}  y{\left (x \right )}$$</div>

<p>To solve with $y(0) = 1$ and $y'(0) = 1$ we have:</p><pre class="sourceCode julia">out = dsolve(eqn, x, (y, 0, 1), (y', 0, 1))</pre>
<div class="well well-sm">
$$y{\left (x \right )} = \left(4 - 3 e^{- x}\right) e^{- 2 x}$$</div>

<p>To make a plot, we only need the right-hand-side of the answer:</p><pre class="sourceCode julia">plot(rhs(out), -1/3, 2)</pre>
    <div id="3f174419-25bc-4717-b0fb-b408c0e0ade2" style="width:600px;height:400px;"></div>
    <script>
    PLOT = document.getElementById('3f174419-25bc-4717-b0fb-b408c0e0ade2');
    Plotly.plot(PLOT, [{"yaxis":"y","y":[-0.3639093211584317,-0.16594525764970514,0.010533477812474485,0.16736061727754267,0.30622781297210305,0.42869511454495834,0.5362006952520834,0.630069880002471,0.7115235244972477,0.7816857912611007,0.8415913651694356,0.8921921481004488,0.9343634695730548,0.968909847655319,0.996570332030307,1.0180234588749337,1.0338918451308037,1.0447464478136919,1.0511105122101012,1.053463231136109,1.052243135876763,1.0478512379756182,1.0406539396961385,1.0309857297226195,1.0191516795015647,1.0054297545390198,0.9900729539595946,0.9733112906935749,0.9553536237847374,0.9363893534987001,0.9165899891556675,0.8961105989082696,0.8750911500312516,0.8536577476815694,0.8319237795218613,0.8099909730743076,0.7879503721828167,0.7658832385067242,0.7438618835463798,0.7219504363078889,0.7002055513488252,0.678677061606992,0.657408580098512,0.6364380542779777,0.6157982765805702,0.5955173544124829,0.5756191426203495,0.5561236412513976,0.5370473612125775,0.5184036602478722,0.5002030514773531,0.48245348657838333,0.46516061553780463,0.44832802476316597,0.4319574552103026,0.416049002063148,0.4006012973889115,0.38561167708706173,0.3710763333533552,0.3569904537899212,0.34334834820865795,0.33014356409745554,0.3173689916466232,0.30501695916595106,0.2930793196607272,0.2815475292774089,0.27041271827619473,0.25966575513817375,0.2492973043687496,0.23929787851641454,0.22965788488642339,0.22036766739228425,0.2114175439540295,0.202797839820762,0.19449891716582107,0.18651120127590198,0.17882520363044896,0.17143154214447434,0.1643209588265077,0.15748433508352255,0.15091270488631048,0.14459726599176698,0.13852938940282122,0.132700627232189,0.12710271912267324,0.12172759736429374,0.1165673908370318,0.11161442789734774,0.10686123831681593,0.102300554372156,0.09792531117757024,0.0937286463425726,0.08970389903136741,0.08584460849326193,0.08214451212753507,0.07859754314059773,0.07519782784813227,0.07193968267015735,0.06881761086259963,0.06582629902493764],"showlegend":true,"name":"y1","type":"scatter","xaxis":"x","line":{"width":1,"dash":"solid","color":"rgba(0, 154, 250, 1.000)","shape":"linear"},"x":[-0.3333333333333333,-0.30976430976430974,-0.28619528619528617,-0.26262626262626265,-0.23905723905723905,-0.21548821548821548,-0.1919191919191919,-0.16835016835016836,-0.1447811447811448,-0.12121212121212122,-0.09764309764309764,-0.07407407407407407,-0.050505050505050504,-0.026936026936026935,-0.003367003367003367,0.020202020202020204,0.04377104377104377,0.06734006734006734,0.09090909090909091,0.11447811447811448,0.13804713804713806,0.16161616161616163,0.18518518518518517,0.20875420875420875,0.23232323232323232,0.2558922558922559,0.27946127946127947,0.30303030303030304,0.3265993265993266,0.3501683501683502,0.37373737373737376,0.39730639730639733,0.4208754208754209,0.4444444444444444,0.468013468013468,0.49158249158249157,0.5151515151515151,0.5387205387205387,0.5622895622895623,0.5858585858585859,0.6094276094276094,0.632996632996633,0.6565656565656566,0.6801346801346801,0.7037037037037037,0.7272727272727273,0.7508417508417509,0.7744107744107744,0.797979797979798,0.8215488215488216,0.8451178451178452,0.8686868686868687,0.8922558922558923,0.9158249158249159,0.9393939393939394,0.9629629629629629,0.9865319865319865,1.0101010101010102,1.0336700336700337,1.0572390572390573,1.0808080808080809,1.1043771043771045,1.127946127946128,1.1515151515151516,1.1750841750841752,1.1986531986531987,1.2222222222222223,1.245791245791246,1.2693602693602695,1.292929292929293,1.3164983164983164,1.34006734006734,1.3636363636363635,1.387205387205387,1.4107744107744107,1.4343434343434343,1.4579124579124578,1.4814814814814814,1.505050505050505,1.5286195286195285,1.5521885521885521,1.5757575757575757,1.5993265993265993,1.6228956228956228,1.6464646464646464,1.67003367003367,1.6936026936026936,1.7171717171717171,1.7407407407407407,1.7643097643097643,1.7878787878787878,1.8114478114478114,1.835016835016835,1.8585858585858586,1.8821548821548821,1.9057239057239057,1.9292929292929293,1.9528619528619529,1.9764309764309764,2.0],"mode":"lines"}], {"yaxis":{"type":"-","titlefont":{"size":15,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"title":"","tickfont":{"size":11,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"ticks":"inside","tickmode":"array","showgrid":true,"tickvals":[0.0,0.5,1.0],"domain":[0.03762029746281716,0.9901574803149606],"ticktext":["0.0","0.5","1.0"],"tickangle":0,"zeroline":false,"linecolor":"rgba(0, 0, 0, 1.000)","tickcolor":"rgba(0, 0, 0, 1.000)","anchor":"x"},"annotations":[],"width":600,"plot_bgcolor":"rgba(255, 255, 255, 1.000)","showlegend":true,"legend":{"bgcolor":"rgba(255, 255, 255, 1.000)","y":1.0,"font":{"size":11,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"bordercolor":"rgba(0, 0, 0, 1.000)","x":1.0},"xaxis":{"type":"-","titlefont":{"size":15,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"title":"","tickfont":{"size":11,"color":"rgba(0, 0, 0, 1.000)","family":"Helvetica"},"ticks":"inside","tickmode":"array","showgrid":true,"tickvals":[0.0,0.5,1.0,1.5,2.0],"domain":[0.05100612423447069,0.9934383202099737],"ticktext":["0.0","0.5","1.0","1.5","2.0"],"tickangle":0,"zeroline":false,"linecolor":"rgba(0, 0, 0, 1.000)","tickcolor":"rgba(0, 0, 0, 1.000)","anchor":"y"},"paper_bgcolor":"rgba(255, 255, 255, 1.000)","height":400,"margin":{"r":0,"l":0,"b":0,"t":20}});
    </script>
<h5>Example</h5><p>Boundary value problems can be solved for as well through a similar syntax. Continuing with examples from the <a href="https://reference.wolfram.com/language/tutorial/DSolveLinearBVPs.html">Wolfram</a> page, we solve $y''(x) +y(x) = e^x$ over $[0,1]$ with conditions $y(0)=1$, $y(1) = 1/2$:</p><pre class="sourceCode julia">eqn = y''(x) + y(x) - exp(x)
dsolve(eqn, x, (y, 0, 1), (y, 1, 1//2))</pre>
<div class="well well-sm">
$$y{\left (x \right )} = \frac{e^{x}}{2} + \frac{\sin{\left (x \right )}}{2 \sin{\left (1 \right )}} \left(- e - \cos{\left (1 \right )} + 1\right) + \frac{1}{2} \cos{\left (x \right )}$$</div>


  </div>
</div>  

</body>
</html>
